Eugene Johnson - Lab 10


% ode_solver.m
% This function first calls RK4 method to get first three
% values of t, this is a jump start because the way the AB/AM
% methods work.

function XX = ode_solver(Xinit,h,f,tfinal);

XX = zeros(5,(tfinal/h)+1);

XX(:,1) = Xinit;

% getting next three from rk4 method 

for i = 1:3,

   XX(:,i+1) = rk4(XX(:,i),h,f);


% Calling ab_am function to get predictor/corrector.

for i = 4:(tfinal/h),

    XX(:,i+1) = ab_am(XX(:,(i-3):i),h,f);



function M = make_movie(XX)     

% make_movie takes as input the XX matrix which is the output of
% ode_solver.  This function should pick about twenty equally-spaced
% columns of XX and make a movie of them.

% clear the figure window and initialize M to be big enough to hold
% fifty frames.  If MATLAB whines about running out of memory at this
% point, make your figure window smaller.

% You really don't get the true picture because the laws of physics
% the physics law is Kepler's law and this causes the shuttle to
% zip around the earth so fast that it looks like it passes through 
% the earth.  I tryed to make it look more reallistic. I got to where
% I wass take 200 pictures and it didn't change the movie so I
% decided that it was best to write a little comment and leave it 
% at 50 clips.  So here we go.

clf reset; M = moviein(50);

% xpos and ypos should end up being 1-by-50 (or smaller) row vectors 
% picked at equal intervals from the second and third row of XX

% ------ finish these lines ------ %

xpos = XX(2,1:(6200/50):6201);
ypos = XX(3,1:(6200/50):6201); 

% determine how to fix the axis limits so the movie doesn't wiggle around
% as it plays

% loop through the values of xpos and ypos.  Plot the earth as a white 'o',
% the moon as a white 'x', and the shuttle as a white '+'.

u = .0121285;

for i = 1:length(xpos),


        axis([min(xpos)-.5 max(xpos)+.5 min(ypos) max(ypos)]); axis equal

        M(:,i) = getframe;




% rk4.m
% This function is used to jump start everything.  It gets the first 
% three values of t

function Z = rk4(X,h,f);

f1 = func(X,f);
f2 = func(X+h/2*f1,f);
f3 = func(X+h/2*f2,f);
f4 = func(X+h*f3,f);

% getting the values of Z for the whole thing to meet

Z = X + h/6*(f1+2*f2+2*f3+f4);


% ab_am.m
% This function uses the AB/AM methods (predictor/corrector) to
% evaluate the values of the differential equations.

% implementing the predictor/corrector methods

function Z = ab_am(XX,h,f);

% Temp is a temperary variable which holds the predicted value
% Getting predicted value.

TEMP = XX(:,4) + h/24*(55*func(XX(:,4),f) - 59*func(XX(:,3),f) + ... 
       37*func(XX(:,2),f) - 9*func(XX(:,1),f));

% Correcting the predicted value

Z = XX(:,4) + h/24*(9*func(TEMP,f) + 19*func(XX(:,4),f) - ... 
    5*func(XX(:,3),f) + func(XX(:,2),f));