Увеличение числа p в соответствии с полученным другим значением h

################################подключение необходимых библиотек####################
import math 
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
##############################определение функции,вычисляющую правую часть системы###
def sir_model(t, y, B, g, R, k):
    v, h, l, teta, p = y
    dv_dt = ((-0.5*B*p*v*v)/2)-((g*math.sin(teta)))
    dh_dt = v*math.sin(teta)
    dl_dt = v*math.cos(teta)
    dteta_dt = ((B*k*v*p)/2)-((g*math.cos(teta))/v)
    dp_dt = p+0.0551
    return dv_dt, dh_dt, dl_dt, dteta_dt, dp_dt
################################исходные данные######################################   
k = .1
B = .006
p0 = .0136
g = 9.81
R = 6371000
v0 = 3000
l0 = 0
h0= 30000
teta0 = 20
###########################решатель##################################################
T = 30
sol = integrate.solve_ivp(lambda t, y: sir_model(t, y, B ,g, R, k), [0,T], [v0,h0,l0,teta0,p0], dense_output = True)
###########################ВИЗУАЛИЗАЦИЯ##############################################
t = np.linspace(0, T, 50)
y = sol.sol(t)
v, h, l, teta, p = y[0,:], y[1,:], y[2,:], y[3,:], y[4,:]
pass