
function [dy] = fun(y)

y_1 = y(1);
y_2 = y(2);
y_3 = y(3);
y_4 = y(4);
y_5 = y(5);
y_6 = y(6);


mu = 1/81.45;

con1 = ((y_2^2)+(y_3^2)+((y_1+mu-1)^2))^(3/2);
con2 = ((y_2^2)+(y_3^2)+((y_1+mu)^2))^(3/2);

dy_1 = y_4;
dy_2 = y_5;
dy_3 = y_6;
dy_4 = (2*y_5)+y_1-((mu*(y_1+mu-1))/con1)-(((1-mu)*(y_1+mu))/con2);
dy_5 = (-2*y_4)+y_2-((mu*y_2)/con1)-(((1-mu)*y_2)/con2);
dy_6 = -((mu*y_3)/con1)-(((1-mu)*y_3)/con2);

dy = [dy_1,dy_2,dy_3,dy_4,dy_5,dy_6];