ode_solver.m

function XX = ode_solver(X, h, f,tfinal)

u = 0.0121285;

XX = zeros(5,(tfinal/h)+1);
XX(:,1) = X;
for n = 2:4,
        XX(:,n) = rk4(XX(:,(n-1)),h,f);
end

for n = 5:(tfinal/h+1),
        XX(:,n) = ab_am(XX(:,[(n-4):(n-1)]),h,f);
end

%min_dist = min(sqrt((XX(2,:)+u).^2+XX(3,:).^2))*238000-4000
hold off; plot(-u,0,'o'); hold on; plot(1-u,0,'x'); plot(XX(2,:),XX(3,:))
axis([-2 2 -2 2]); axis equal

make_movie.m

function M = make_movie(XX)     

% make_movie takes as input the XX matrix which is the output of
% ode_solver.  This function should pick about twenty equally-spaced
% columns of XX and make a movie of them.

% clear the figure window and initialize M to be big enough to hold
% fifty frames.  If MATLAB whines about running out of memory at this
% point, make your figure window smaller.

clf reset; M = moviein(50);

% xpos and ypos should end up being 1-by-50 (or smaller) row vectors 
% picked at equal intervals from the second and third row of XX

% ------ finish these lines ------ %

[m,n]=size(XX);
xpos = XX(2,1:((n-1)/50):n);
ypos = XX(3,1:((n-1)/50):n);

% determine how to fix the axis limits so the movie doesn't wiggle around
% as it plays

% loop through the values of xpos and ypos.  Plot the earth as a white 'o',
% the moon as a white 'x', and the shuttle as a white '+'.

u = .0121285;

for i = 1:length(xpos),

        hold off; plot(-u, 0, 'o'); hold on; plot(1-u, 0, 'x'); plot(xpos(i), ypos(i), '+');
axis([-2 2 -2 2]); axis equal
%axis([min(xpos) max(xpos) min(ypos) max(ypos)]); axis equal
        M(:,i) = getframe;


end

rk4.m

function Z = rk4(X, h, f)

F1 = func(X,f);
F2 = func(X + .5*h*F1,f);
F3 = func(X + .5*h*F2,f);
F4 = func(X + h*F3,f);

Z = X + h/6*(F1 + 2*F2 + 2*F3 + F4);

ab_am.m

Z = XX(:,4) + h/24 * (55 * func(XX(:,4),f) - 59* func(XX(:,3),f) + 37 * func(XX(:,2),f) - 9*func(XX(:,1),f));

Z = XX(:,4) + h/24 * (9*func(Z, f) + 19 * func(XX(:,4),f) - 5* func(XX(:,3),f) + func(XX(:,2),h));