function varargout=DYDt3(t,Y,flag,G,ME,tspan,Y0)
% DYDt3.m 	供主程序調用的ODE函數文件
%					本文件自帶三個子函數：f,fi,fev。
%		t, Y		分別是自變量和一階函數向量，是最基本的輸入宗量。
%		flag		第三輸入宗量，它專供解算指令（如ode45）作調用通知。
%					在運行中，解算指令會根據需要向flag發不同的字符串。
% varargout 是"變"輸出宗量。它由變維的元胞數組構成。每個元胞中可以存放指令
%					所產生的任意形式的數據。
switch flag
case ''      	% 必須用空串符。情況為"真"時，計算導數 dY/dt = f(t,Y)。
  varargout{1} = f(t,Y,G,ME);		%輸出為一個元胞，容納f子函數的一個輸出Yd
case 'init'  	% 必須用'init'，情況為"真"時，傳送計算區間、初值、設置參數。
  [varargout{1:3}] = fi(tspan,Y0);	%輸出為三個元胞，容納fi子函數的三個輸出量
case 'events'	% 必須用'events'，情況為"真"時，設置事件性質。
  [varargout{1:3}] = fev(t,Y,Y0);	%輸出三個元胞，容納fev子函數的三個輸出量
otherwise
  error(['Unknown flag ''' flag '''.']);
end
% ------------------------------------------------------------------
function Yd = f(t,Y,G,ME)
%計算導數子函數，被"父"函數DYDt3調用。
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)
%設置時間區間、初值、算法參數子函數，被"父"函數DYDt3調用。
ts=tspan;y0 = Y0;
options = odeset('Events','on','Reltol',1e-5,'Abstol',1e-4);
						%開動ode45的"事件判斷"功能，設置相對誤差和絕對誤差。
% ------------------------------------------------------------------
function [value,isterminal,direction] = fev(t,Y,Y0)
%事件判斷子函數，被"父"函數DYDt3調用。
dDSQdt = 2 * ((Y(1:2)-Y0(1:2))' * Y(3:4));
%dDSQdt是離初始點的二乘距離u的時間導數du/dt，而u=(x-x0)^2+(y-y0)^2 。
value = [dDSQdt; dDSQdt];	%定義兩個穿越0的事件
direction = [1; -1];	%第一事件：以漸增方式穿越0。第二事件：以漸減方式穿越0
isterminal = [1; 0];	%第一事件發生后，終止計算；而第二事件發生后，繼續計算。