% 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); end % Calling ab_am function to get predictor/corrector. for i = 4:(tfinal/h), XX(:,i+1) = ab_am(XX(:,(i-3):i),h,f); end

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), plot(-u,0,'wo',1-u,0,'wx',xpos(i),ypos(i),'w+'); axis([min(xpos)-.5 max(xpos)+.5 min(ypos) max(ypos)]); axis equal M(:,i) = getframe; end movie(M,2);

% 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));