ode45求解微分方程

%%% 主函数
clc
clear
close all


y(1)=5600;
y(2)=-pi/12;
y(3)=6571000;
y(4)=0;

tspan=[0:1:160];
y0=y';
[t,res]=ode45(@myode1,tspan,y0);
x11= res(:,1);
x12=res(:,2);
x13=res(:,3);
x14=res(:,4);
dv=diff(v);
%% plot
% Draw;			%绘图函数，自己定义即可

ode45函数

function dY = myode1(t, Y)
dY=zeros(4,1);
R_Earth=6371*1000;             
pp = 0.2;
S_M= 1/4*pi*(1.67^2); 

m = 955;
v = Y(1);
theta = Y(2);

r = Y(3);
beta_e = Y(4);

h = r-R_Earth;                  
g=9.8;
rho=1.22;
q = 0.5*rho*v^2;
if (h>=0)
dY(1)=-pp*q*S_M/m - g * sin(theta);
dY(2)= (v /r - g / v) * cos(theta);
dY(3)=v*sin(theta);
dY(4)=v/r*(cos(theta));
end

end