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

a = 1; b = 1; c = 0;

% Part (a) one set of initial conditions: 
%
x_0 = [0.1;0]; % the initial condition
t_0 = 0.0; % the initial time
nSteps = 20;
t_final = 2.0;
tau = t_final/nSteps; % the timestep 
params = [ a, b, c ];

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, @van_der_pol_eq, params );
  t_history(ii+1) = tim1+tau;
  x_history(ii+1,:) = xi(:)';
end
x_history_pt_a = x_history;

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


% Part (b) another set of initial conditions (code from the above copied): 
%
x_0 = [0.2;0]; % the initial condition

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, @van_der_pol_eq, params );
  t_history(ii+1) = tim1+tau;
  x_history(ii+1,:) = xi(:)';
end
x_history_pt_b = x_history;

figure(fh); 
ph_b = plot( t_history, x_history(:,1), '-x' ); 


% Part (d) integrate the LINEARIZED van-der-pol equation:
%
x_0 = [0.1;0]; % the initial condition

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, @van_der_pol_eq_linearized, params );
  t_history(ii+1) = tim1+tau;
  x_history(ii+1,:) = xi(:)';
end
x_history_pt_d = x_history;

figure(fh); 
ph_d = plot( t_history, x_history_pt_a(:,1)+x_history_pt_d(:,1), 'r-x' ); % add the linearized solution to that from part_a


xlabel('time'); ylabel('x(t) for the van der pol equation');
axis('tight');

legend([ph_a,ph_b,ph_d], {'NL: x(0)=0.1','NL: x(0)=0.2','L: x(0)+\Delta x(0)'}, 'location', 'northeast');

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