function [yp]=f1(t,y)
yp=[];
yp = -y;


function [yp]=f2(t,y)
yp=[];
yp = -50*(y-cos(t));


function [ydot] = aren(t, y)

mu = 0.012277471;
r1 = sqrt((y(1)+mu)^2 + y(2)^2);
r2 = sqrt((y(1)-1+mu)^2 + y(2)^2);

ydot(1) = y(3);
ydot(2) = y(4);
ydot(3) = y(1) + 2*y(4) - (1-mu)*(y(1)+mu)/r1^3 - mu*(y(1)-1+mu)/r2^3;
ydot(4) = y(2) - 2*y(3) - (1-mu)*y(2)/r1^3 - mu*y(2)/r2^3;