function [du] = Apollosys(t,u)
mu=1/82.45;          % Ratio of masses: M_moon/M_earth;  
lam=1-mu;
% Auxiliary variables:
% x=u1,x'=u2,y=u3,y'=u4
% For convenience take local variables:
u1=u(1);u2=u(2);u3=u(3);u4=u(4);
r1=sqrt((u1+mu).^2 + u3.^2);
r2=sqrt((u1-lam).^2 + u3.^2);
% Note: all variables are scalars =>  no need for .^ (works though).
du=[u2;
    2*u4+u1-lam*(u1+mu)./r1.^3-mu*(u1-lam)./r2.^3;
    u4;
    -2*u2+u3-lam*u3./r1.^3-mu*u3./r2.^3];
end

