Отображение графика в решении задачи с модулем SciPy

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

pass

plt.show()