用python求解二阶常系数微分方程数值解

import matplotlib.pyplot as plt
from scipy import linspace,exp
from scipy.integrate import odeint
import numpy as np
import mpl_toolkits.axisartist as axisartist
import matplotlib
#创建画布
fig=plt.figure()
#添加画布
ax=axisartist.Subplot(fig,111)
fig.add_axes(ax)
#去掉右上角的边界线
ax.axis['top'].set_visible(False)
ax.axis['right'].set_visible(False)
#给坐标轴设置箭头和大小
ax.axis['bottom'].set_axisline_style("->",size=1.5)
ax.axis['left'].set_axisline_style("->",size=1.5)
#设置字体样式和风格
matplotlib.pyplot.rcParams['font.family']='serif'
matplotlib.pyplot.rcParams['font.serif']=['Times New Roman']+matplotlib.pyplot.rcParams['font.serif']
matplotlib.pyplot.rcParams['font.size']=14
def fun(t,y):
    dy1=y[1]#一阶导数
    dy2=9.8-(1/50)*0.01*y[1]**2#二阶导数
    return [dy1,dy2]#返回一阶导和二阶导的表达式
def solve():
    y0=[300,0]#300表示题中给定的一阶导初值，0表示二阶导初值
    t=linspace(0,20,100)#给出自变量的取值范围
    y=odeint(fun,y0,t,tfirst=True)
    y1, = ax.plot(t, y[:, 0], label='y')#画出y的值对应的函数图像
    y1_1, = ax.plot(t, y[:, 1], label="y'")#画出y的一阶导对应的函数图像
    plt.xlabel("t")
    plt.ylabel("y")#设置标签
    plt.legend(handles=[y1, y1_1])
    plt.show()
    print(y1_1[:,0])#原函数的数值解
    print(y1[:,1])#一阶导的数值解
solve()