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