% right-hand side of arenstorf problem function yDot = arenstorf(t,y) global mu muHat % unpack y u1 = y(1); u2 = y(2); u1Dot = y(3); u2Dot = y(4); % define variables of convenience r1 = ((u1+mu)^2+u2^2)^(3/2); r2 = ((u1-muHat)^2+u2^2)^(3/2); u1DDot = u1 + 2*u2Dot - muHat*(u1+mu)/r1 - mu*(u1-muHat)/r2; u2DDot = u2 - 2*u1Dot - muHat*u2/r1 - mu*u2/r2; % pack up derivatives yDot = [u1Dot u2Dot u1DDot u2DDot]';