Matlab Euler's method

% matlab script to test efficiency of
% Euler's method, classical Runge-Kutta, and ode45
% on Arenstorf orbit problem

close all
clear all

% these are variables we would like the right-hand function to "see"
% without actually passing them as arguments
global mu muHat

% set normalized masses
mu = 0.012277471;
muHat = 1 - mu;

% set time span of integration
t0 = 0;
tf = 17.1;

% set initial conditions
u1 = 0.994;
u2 = 0;
u1Dot = 0;
u2Dot = -2.00158510637908252240537862224;

% pack into vector y0
y0 = [u1 u2 u1Dot u2Dot]';

% set the name of the function to compute the right-hand side
f = 'arenstorf';
disp(['The ' f ' problem'])
disp(' ')

% begin with Euler's method
disp('Experiments with Euler''s method.')

% "logical" variable to end given phase (Euler, Runge-Kutta)
iMore = 1;

while (iMore == 1),
    % prompt user for number of steps
    nSteps = input('Please enter the number of steps to take: ');
    % set up array to store solution
    y = zeros(length(y0),nSteps+1);
    h = (tf-t0)/nSteps;
    t = linspace(t0,tf,nSteps+1);
    y(:,1) = y0;
    tic
    for i=1:nSteps,
        y(:,i+1) = y(:,i) + h*feval(f,t(i),y(:,i));
    end
    toc
    plot(y(1,:),y(2,:))
    iMore = input('Do you wish to repeat with Euler''s method? (1=yes,0=no) ');
    disp(' ')
end
    
close all
% now investigate classical Runge-Kutta 
disp('Experiments with classical Runge-Kutta method.')

iMore = 1;

while (iMore == 1),
    % prompt user for number of steps
    nSteps = input('Please enter the number of steps to take: ');
    % set up array to store solution
    y = zeros(length(y0),nSteps+1);
    h = (tf-t0)/nSteps;
    t = linspace(t0,tf,nSteps+1);
    y(:,1) = y0;
    tic
    for i=1:nSteps,
        K1 = feval(f,t(i),y(:,i));
        K2 = feval(f,t(i)+h/2,y(:,i)+h*K1/2);
        K3 = feval(f,t(i)+h/2,y(:,i)+h*K2/2);
        K4 = feval(f,t(i+1),y(:,i)+h*K3);
        y(:,i+1) = y(:,i) + h*(K1+2*(K2+K3)+K4)/6;
    end
    toc
    plot(y(1,:),y(2,:))
    iMore = input('Do you wish to repeat with classical RK? (1=yes,0=no) ');
    disp(' ')
end
close all


% finally use ode45
disp('Solving problem with ode45 ... ')

% try with and without options in call to ode45
options = odeset('RelTol',1e-4);

tic
[t,y] = ode45(f,[t0 tf],y0,options);
toc
plot(y(:,1),y(:,2))

% from Wikipedia about the Arenstorf orbit

% Richard F. Arenstorf is an American mathematician who discovered a
% stable orbit between the Earth and the Moon, called an Arenstorf
% Orbit, which was the basis of the orbit used by the Apollo Program
% for going to the Moon.

% The Arenstorf Orbit

% While the orbit of a satellite around the Sun was empirically
% discovered by Kepler and theoretically proven by Newton to be an
% ellipse, at the time when the United States was interested in going
% to the Moon, there was no such solution known for the shape of a
% satellite orbiting regularly around two objects, such as a
% spacecraft going between the Earth and the Moon. This is a special
% case of the infamous Three Body Problem, for which a general
% analytical solution is not known because of its complexity of
% solving the effect of three bodies which all pull on each other
% while moving, a total of six interactions. However the case of an
% Earth-Moon satellite can be simplified to four interactions, because
% although the three objects gravitationally all pull on each other,
% the effect of the spacecraft's gravity upon the motion of the vastly
% more massive Earth and Moon is practically non-existent. Arenstorf
% found a stable orbit for a spacecraft orbiting between the Earth and
% Moon, shaped like an '8' with the Earth or Moon located inside each
% loop of the '8'. This orbit is the basis of a path going to the Moon
% from the Earth, such as the United States Apollo program. For a
% permanent presence on the Moon, it would be the path of what
% Arenstorf calls a 'Space Bus', a ferry which could regularly orbit
% supplies and people between the Earth and Moon without directly
% expending fuel. By staying on the Arenstorf orbit, lunar astronauts
% automatically return back to Earth. Before leaving NASA at the first
% Moon Landing, Arenstorf mapped out an emergency rescue orbit, which
% was used in the Apollo 13 incident, in which a catastrophic
% malfunction forced aborting the Moon landing, but the astronauts
% ultimately returned safely to Earth without a major course
% adjustment.