X0 = [0; 0;         % Sun pos
     1.496e11; 0;   % Earth pos
     1.496e11+3.84403e8; 0;    % Moon pos
     0; 0;          % Sun vel
     0; 2.98e4;     % Earth vel
     0; 2.98e4+1.03e3      % Moon vel
];

% Step through one year
year = 365*24*60*60;
options = odeset('RelTol',1e-5);
[T, X] = ode45(@f, [0 year], X0, options);

% Plot the earth's track
clf
plot(X(:, 3), X(:,4), 'b')
hold on;
axis square;

% Plot the moon's track
plot(X(:, 5), X(:,6), 'g')

% Mark the origin
plot(0, 0, 'r*');


[m, n] = size(X)

for i=1:m
    clf

    earthPosX = X(i,3);
    earthPosY = X(i,4);
 
    plot(X(1:i, 3), X(1:i,4), 'b')
    hold on;
    axis square;
    plot(X(i, 3), X(i,4), 'b*')
    plot(X(1:i, 5), X(1:i,6), 'g')
    plot(X(i, 5), X(i,6), 'g*')

    axis([earthPosX-1e9, earthPosX+1e9 earthPosY-1e9, earthPosY+1e9]);
    A(i) = getframe;
end
