clear;
y0=0;
x0=0;                 						%初始位置
vMag=input('輸入初始速度(m/s)：');        		%輸入初始速度
vDir=input('輸入初始方向(度)');
tf=input('輸入飛行時(shí)間(秒)');            			%輸入飛行時(shí)間
vx0=vMag*cos(vDir*(pi/180));           			%計(jì)算x，y方向的初始速度
vy0=vMag*sin(vDir*(pi/180));
r0=[0;0;vx0;vy0];
[t,r]=ode45('e31a',[0,tf],r0),           			%數(shù)值積分（調(diào)用函數(shù)程序e31a.m）
plot(r(:,1),r(:,2)),hold on                   		%計(jì)算軌跡
%ode45規(guī)定返回的結(jié)果中：t是列向量，各時(shí)刻的r為4列向量
%注意下一個(gè)語(yǔ)句的意義：找y<0的下標(biāo)所對(duì)應(yīng)的x的最小值，以粗略計(jì)算射程
xmax=min(r(find(r(:,2)<0),1))
plot([0,150],[0,0])                    			%畫(huà)出x坐標(biāo)線