global G ME						%在主程序中定義全局變量傳遞參數		<1>
G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9;
tspan=[t0,tf];					%指定解算微分方程的時間區(qū)間
Y0=[x0;0;0;vy0];					%按式(4.14.2.1-6)給定初值向量
[t,YY]=ode45('DYDt',tspan,Y0);										%<8>
X=YY(:,1);						%輸出Y的第一列是位移數據x(t)
Y=YY(:,2);						%輸出Y的第二列是位移數據y(t)
plot(X,Y,'b','Linewidth',2); hold on
axis('image')						%保證x、y軸等長刻度，且坐標框恰包容圖形
[XE,YE,ZE] = sphere(10);		%產生單位球面數據
RE=0.64e7;						%地球半徑
XE=RE*XE;YE=RE*YE;ZE=0*ZE;		%坐標紙上的地球平面數據
mesh(XE,YE,ZE),hold off			%繪地球示意圖