ошибка при решении системы дифференциальных уравнений с помощью SciPy.PYTHON

необходимо решить систему дифференциальных уравнений и построить график с помощью модулей numpy, matplotlib.pyplot, scipy. Есть система уравнений рисунок 1 и входные данные рисунок 2. Ввел вё как надо но выдает ошибку.

ОШИБКА:

Traceback (most recent call last):
  File "C:\Users\srvadm\Desktop\еще.py", line 31, in <module>
    plt.plot(t,v)
  File "C:\Users\srvadm\AppData\Local\Programs\Python\Python310\lib\site-packages\matplotlib\pyplot.py", line 2740, in plot
    return gca().plot(
  File "C:\Users\srvadm\AppData\Local\Programs\Python\Python310\lib\site-packages\matplotlib\axes\_axes.py", line 1662, in plot
    lines = [*self._get_lines(*args, data=data, **kwargs)]
  File "C:\Users\srvadm\AppData\Local\Programs\Python\Python310\lib\site-packages\matplotlib\axes\_base.py", line 311, in __call__
    yield from self._plot_args(
  File "C:\Users\srvadm\AppData\Local\Programs\Python\Python310\lib\site-packages\matplotlib\axes\_base.py", line 504, in _plot_args
    raise ValueError(f"x and y must have same first dimension, but "
ValueError: x and y must have same first dimension, but have shapes (50,) and (4, 50)

МОЙ КОД:

################################подключение необходимых библиотек####################
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 = ((-B*v*v*p)/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.sin(teta))/v)*((1-v*v)/(g*R))
    return dv_dt, dh_dt, dl_dt, dteta_dt
################################исходные данные######################################   
k = 3.2
B = .006
p = .225
g = 9.81
R = 6371000
v0 = 7000
l0 = 0
h0= 80000
teta0 = -0.349
###########################решатель##################################################
T = 50
sol = integrate.solve_ivp(sir_model, [0,T], [v0,h0,l0,teta0], dense_output = True,args = [k,B,p,g,R])
###########################ВИЗУАЛИЗАЦИЯ##############################################
t = np.linspace(0, T)
y = sol.sol(t)
v, h, l, teta = y[0 :], y[1, :], y[2, :], y[3, :]
plt.plot(t,v)
plt.plot(t,h)
plt.plot(t,l)
plt.plot(t,teta)
pass

Подскажите, пожалуйста, как мне получить график