    function z = fd(n,x,y)

    mu = 0.012277471;
    muhat = 1 - mu; 
    D1 = ((y(1) + mu)^2 + y(2)^2)^1.5;
    D2 = ((y(1) - muhat)^2 + y(2)^2)^1.5;
    z = [y(3); y(4); y(1)+2*y(4)-muhat*(y(1)+mu)/D1-mu*(y(1)-muhat)/D2; y(2)-2*y(3)-muhat*y(2)/D1-mu*y(2)/D2];
