%   
% Written by:
% -- 
% John L. Weatherwax                2006-08-28
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; clc; clear; 

% for RK4:
addpath('../../../garcia/edition2/Matlab'); 


% Part (a) the integraion based on the function for of f() in dxdt = f(x,t) :
%
x_0 = [0;0]; % the initial conditions
t_0 = 0.0; % the initial time
nSteps = 100;
t_final = 10;
tau = t_final/nSteps; % the timestep 

Phi_DT = [ exp(-0.5*tau), 5*(exp(-0.5*tau) - exp(-0.7*tau)); 0, exp(-0.7*tau) ];
Gamma_DT = -(10/7)*[ 7*exp(-0.5*tau) - 5*exp(-0.7*tau) - 2; exp(-0.7*tau)-1 ];



t_history = zeros(nSteps+1,1);
x_history = zeros(nSteps+1,2);
t_history(1) = t_0; 
x_history(1,:) = x_0(:).';
for ii=1:nSteps,
  xim1 = x_history(ii,:)';
  tim1 = t_history(ii);
  
  xi = Phi_DT * xim1 + Gamma_DT * 1; % the integration step 
  
  t_history(ii+1) = tim1+tau;
  x_history(ii+1,:) = xi(:)';
end


fh = figure; 
ph_lin = plot( t_history, x_history(:,1), '-xr' ); hold on; 


% Part (b) the RK4 integrations:
%
params = [];

t_history = zeros(nSteps+1,1);
x_history = zeros(nSteps+1,2);
t_history(1) = t_0; 
x_history(1,:) = x_0(:).';
for ii=1:nSteps,
  xim1 = x_history(ii,:)';
  tim1 = t_history(ii);
  xi   = rk4( xim1, tim1, tau, @sect_2_3_prob_4_eq, params );
  t_history(ii+1) = tim1+tau;
  x_history(ii+1,:) = xi(:)';
end

figure(fh);
ph_rk4 = plot( t_history, x_history(:,1), '-og' ); hold on; 

legend( [ph_lin,ph_rk4], {'special linear integrator', 'RK4'}, 'location', 'northwest' );
xlabel('time'); ylabel('x(t)');

saveas(gcf,'../../WriteUp/Graphics/Chapter2/sect_2_3_prob_4.eps','epsc');