function dydt = fct_tbp(t,y)
% Refer to 5-48,5-49,5-50, "Matlab, Mathematic. The language of Technical 
% Computing". The MathWorks. 2004. Matlab7. 
% The system is defined as follows:
% (1) y_1'=y_3
% (2) y_2'=y_4
% (3) y_3'=2(y_4)+(y_1)-[(µ*((y_1)+µ)/(r_1)^3]-[(µ(y_1)+µ*)/(r_2)^3]
% (4) y_4'=-2(y_3)+(y_2)-[µ*(y_2)/(r_1)^3]-[(µ(y_2)/(r_2)^3]
% with 
% µ=1/82.45
% µ*=1-µ
% r1=[(y_1)+µ)^2+(y_2)^2]^(1/2) 
% r2=[(y_1)-µ*)^2+(y_2)^2]^(1/2)

mu=1/82.45;
mustar=1-mu;
r13=((y(1)+mu)^2+y(2)^2)^1.5;
r23=((y(1)-mustar)^2+y(2)^2)^1.5;
dydt= [ y(3)
        y(4)
        2*y(4)+ y(1)- mustar*((y(1)+mu)/r13)- ...
                      mu*((y(1)-mustar)/r23)
        -2*y(3)+y(2)-mustar*(y(2)/r13)-mu*(y(2)/r23) ];
        
% The MathWorks, Matlab7. Copyright (c).