function t=orbitode113jac

% Problem parameters
mu = 1 / 82.45; % Shampine & Gordon
mustar = 1 - mu;
y0 = [1.2; 0; 0; -1.04935750983031990726; 0; 0; 0; 0]; % Shampine & Gordon
tspan = [0 7]; % Shampine & Gordon
options = odeset('RelTol',1e-8,'AbsTol',1e-10,'Events',@events);

[t,y,te,ye,ie] = ode113(@f,tspan,y0,options);

% Plotting circles representing earth and moon
rade = 6371/384403;
radm = 1737/384403;
th=linspace(-pi,pi,101);
esx = rade*cos(th)-mu;
esy = rade*sin(th);
msx = radm*cos(th)+1-mu;
msy = radm*sin(th);
close all
figure
% Plot timesteps scaled by speed so they are "distance-steps"
avevel = (y(1:end-1,3:4)+ y(2:end,3:4) )/2;
avespeed = sqrt( avevel(:,1).^2 + avevel(:,2).^2 );
%Plots Below: 'k--' line is the distance of the satellite to the moon; 'k' line is the time-steps scaled by average speed; 'k-.' line is the distance of the satellite to the earth
semilogy( t, sqrt( (y(:,1)-mustar).^2 + y(:,2).^2 ), 'k--',t(2:end), diff(t).*avespeed, 'k', t, sqrt( (y(:,1)+mu).^2 + y(:,2).^2 ), 'k-.' ),set(gca,'fontsize',16)
axis([t(1),t(end),1.0e-4,1.0e1]),xlabel('t','fontsize',18),ylabel('distance/stepsize','fontsize',18)
figure
plot(y(:,1),y(:,2),'k-.',esx,esy,'k',msx,msy,'k'),set(gca,'fontsize',18)
hold on
% Estimate of the size of filled circle was obtained by trial and error on actual radius of earth and moon figures, later not needed.
scatter( -mu, 0, 50, 'k', 'filled');
scatter( 1-mu, 0, 50*(radm/rade)^2, 'k', 'filled' );
axis([-1.5,1.5,-1.5,1.5])
axis('square')
xlabel('x','fontsize',18),ylabel('y','fontsize',18)

hold off
figure
semilogy(t(2:end),diff(t), 'k--', t,sqrt(y(:,5).^2+y(:,6).^2+y(:,7).^2+y(:,8).^2),'k'),set(gca,'fontsize',18)
xlabel('t','fontsize',18),ylabel('sensitivity/stepsizes','fontsize',18)
axis([t(1),t(end),1.0e-5,1.0e5])
figure

% Inefficiently, solve it again, so we may compute the residual (should have done this the first time)
sol = ode113(@f,tspan,y0,options);
tt = RefineMesh( sol.x, 10 ); 
np = length(tt);
[yb,ypb] = deval( sol, tt );
res = zeros(1,np);
speed = zeros(size(res));
for i=1:np,
    rr = ypb(:,i) - f(tt(i), yb(:,i) );
    res(i) = norm(rr,inf); speed(i) = norm( ypb(:,i), inf );
end;
semilogy( tt, res./speed, 'k.' ),set(gca,'fontsize',18)
axis([tt(1), tt(end), 1.0e-10,1.0e-7])
xlabel('t','fontsize',18),ylabel('scaled residual','fontsize',18)
  % -----------------------------------------------------------------------
  % Nested functions -- problem parameters provided by the outer function.
  %
  
  function dydt = f(t,y)
  % Derivative function -- mu and mustar shared with the outer function.
    r13 = ((y(1) + mu)^2 + y(2)^2) ^ 1.5;
    r23 = ((y(1) - mustar)^2 + y(2)^2) ^ 1.5;

    % Jacobian computed by Maple
    cg0 = zeros(4,4);
    cg2 = zeros(4,1);
    t1 = 1 - mu;
t2 = y(1) + mu;
t3 = t2 ^ 2;
t4 = y(2) ^ 2;
t5 = t3 + t4;
t6 = sqrt(t5);
t9 = t1 / t6 / t5;
t10 = t1 * t2;
t11 = t5 ^ 2;
t13 = 0.1e1 / t6 / t11;
t14 = 0.2e1 * t13 * t2;
t18 = mu / t1;
t19 = y(1) - 1 + mu;
t20 = t19 ^ 2;
t21 = t20 + t4;
t22 = sqrt(t21);
t25 = t18 / t22 / t21;
t26 = t21 ^ 2;
t28 = 0.1e1 / t22 / t26;
t29 = t19 * t28;
cg0(1,1) = 0;
cg0(1,2) = 0;
cg0(1,3) = 1;
cg0(1,4) = 0;
cg0(2,1) = 0;
cg0(2,2) = 0;
cg0(2,3) = 0;
cg0(2,4) = 1;
cg0(3,1) = (0.1e1 - t9 + 0.3e1 / 0.2e1 * t10 * t14 - t25 + 0.3e1 * t18 * t29 * t19);
cg0(3,2) = (0.3e1 * t10 * t13 * y(2) + 0.3e1 * t18 * t29 * y(2));
cg0(3,3) = 0;
cg0(3,4) = 2;
cg0(4,1) = (0.3e1 / 0.2e1 * t1 * y(2) * t14 + 0.3e1 * t18 * y(2) * t28 * t19);
cg0(4,2) = (0.1e1 - t9 + 0.3e1 * t1 * t4 * t13 - t25 + 0.3e1 * t18 * t4 * t28);
cg0(4,3) = -2;
cg0(4,4) = 0;

% df/dmu
t1 = y(1) + mu;
t2 = t1 ^ 2;
t3 = (y(2) ^ 2);
t4 = t2 + t3;
t5 = sqrt(t4);
t7 = 0.1e1 / t5 / t4;
t9 = 1 - mu;
t12 = t4 ^ 2;
t15 = 0.2e1 / t5 / t12 * t1;
t18 = 1 / t9;
t19 = y(1) - 1 + mu;
t21 = t19 ^ 2;
t22 = t21 + t3;
t23 = sqrt(t22);
t25 = 0.1e1 / t23 / t22;
t27 = t9 ^ 2;
t29 = mu / t27;
t32 = mu * t18;
t34 = t22 ^ 2;
t36 = 0.1e1 / t23 / t34;
cg2(1) = 0;
cg2(2) = 0;
cg2(3) = (t1 * t7 - t9 * t7 + 0.3e1 / 0.2e1 * t9 * t1 * t15 - t18 * t19 * t25 - t29 * t19 * t25 - t32 * t25 + 0.3e1 * t32 * (t19 ^ 2) * t36);
cg2(4) = (y(2) * t7 + 0.3e1 / 0.2e1 * t9 * y(2) * t15 - t18 * y(2) * t25 - t29 * y(2) * t25 + 0.3e1 * t32 * y(2) * t36 * t19);

    cg2 = cg0*y(5:8) + cg2;

    dydt = [ y(3)
             y(4)
             2*y(4) + y(1) - mustar*((y(1)+mu)/r13) - mu/mustar*((y(1)-mustar)/r23)
            -2*y(3) + y(2) - mustar*(y(2)/r13) - mu/mustar*(y(2)/r23)
            cg2(1)
            cg2(2)
            cg2(3)
            cg2(4)];
  end  

  % -----------------------------------------------------------------------

  function [value,isterminal,direction] = events(t,y)
  % Event function -- y0 shared with the outer function.
  % Locate the time when the object returns closest to the initial point y0
  % and starts to move away, and stop integration.  Also locate the time when
  % the object is farthest from the initial point y0 and starts to move closer.
  % 
  % The current distance of the body is
  % 
  %   DSQ = (y(1)-y0(1))^2 + (y(2)-y0(2))^2 = <y(1:2)-y0(1:2),y(1:2)-y0(1:2)>
  %   
  % A local minimum of DSQ occurs when d/dt DSQ crosses zero heading in
  % the positive direction.  We can compute d/dt DSQ as
  % 
  %   d/dt DSQ = 2*(y(1:2)-y0)'*dy(1:2)/dt = 2*(y(1:2)-y0)'*y(3:4)
  %   
  % y0 is shared with the outer function.
  
    dDSQdt = 2 * ((y(1:2)-y0(1:2))' * y(3:4)); 
    value = [dDSQdt; dDSQdt];
    %isterminal = [1;  0];         % stop at local minimum
    isterminal = [0; 0];
    direction  = [1; -1];         % [local minimum, local maximum]
  end

  % -----------------------------------------------------------------------

end  % orbitode

  
  