function varargout=DYDt4(t,Y,flag,G,ME,tspan,Y0)
switch flag
case ''
   varargout{1} = f(t,Y,G,ME);		
case 'init'  
  [varargout{1:3}] = fi(tspan,Y0);
otherwise
  error(['Unknown flag ''' flag '''.']);
end
% ------------------------------------------------------------------
function Yd = f(t,Y,G,ME)
X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.^2));Yd=[V; -G*ME*X/r^3];
% ------------------------------------------------------------------
function [ts,y0,options] = fi(tspan,Y0)
ts=tspan;y0 = Y0;
% 采用向域直接賦值法，設(shè)置options屬性。以供與odeset使用方法對(duì)照。
options.RelTol=1e-5;options.AbsTol=1e-4;
options.OutputFcn='odephas2';	%在積分進(jìn)程中，繪制相平面圖。
options.OutputSel=[1 3];%解向量的第1、3分量分別為相平面圖的橫、縱坐標(biāo)量。