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

close all; clc; clear; 

addpath('../../BookCode');

% Problem EPage 118; Shooting method on EPage 106; 
% 

% Part (i):
%
s = -0.5:0.1:0.5; % evalute the value of y(1) at each of these slopes at zero i.e. y'(0) (we want y(1) ~ 2)
ncases = length(s); b = zeros(1,ncases);
for i=1:ncases,
  [x,y] = ode45( 'c_6_p_3_ode_fn_part_i', [0, 1], [0, s(i)] );
  [m,n] = size(y);
  b(1,i) = y(m,1); % save the value of y(1) for this slope 
end
s0 = aitken(b,s,2); % use aitken's method to find the slope that gives y(1) == 2
[x,y] = ode45( 'c_6_p_3_ode_fn_part_i', [0, 1], [0, s0] );
plot( x, y(:,1), '-k' ); hold('on');

for n=[10,50] 
  x = linspace( 0, 1, n+1 );
  C = ones( 1, n+1 );
  D = -62 * ones( 1, n+1 );
  E = 120 * ones( 1, n+1 ); 
  F = zeros( 1, n+1 ); 
  y_tp = twopoint(x,C,D,E,F,1,1,0,2);

  if n==10
    plot( x, y_tp, ':r' ); 
  else
    plot( x, y_tp, '-r' );
  end
end

xs = linspace( 0, 1, 100 );
ys = ( 1.7513021525390304 * 10^(-26) ) * ( exp( 60 * xs ) - exp( 2 * xs ) ); 
plot( xs, ys, '-g' );
xlabel('x'); ylabel('y(x)'); title('Ex 6.2.i');

legend( 'shooting method', 'twopoint (n=10)', 'twopoint (n=50)', 'exact', 'location', 'best' );

%saveas( gcf, '../../WriteUp/Graphics/Chapter6/c_6_p_3_plot_i.eps', 'epsc' );

% Part (ii):
%
figure;

s = 0:-30:-150; % evalute the value of y(1) at each of these slopes at zero i.e. y'(0) (we want y(1) ~ 0)
ncases = length(s); b = zeros(1,ncases);
for i=1:ncases,
  [x,y] = ode45( 'c_6_p_3_ode_fn_part_ii', [0, 1], [2, s(i)] );
  [m,n] = size(y);
  b(1,i) = y(m,1); % save the value of y(1) for this slope 
end
s0 = aitken(b,s,0); % use aitken's method to find the slope that gives y(1) == 0
[x,y] = ode45( 'c_6_p_3_ode_fn_part_ii', [0, 1], [2, s0] );
plot( x, y(:,1), '-k' ); hold('on');

for n=[10,50] 
  x = linspace( 0, 1, n+1 );
  C = ones( 1, n+1 );
  D = 62 * ones( 1, n+1 );
  E = 120 * ones( 1, n+1 ); 
  F = zeros( 1, n+1 ); 
  y_tp = twopoint(x,C,D,E,F,1,1,2,0);

  if n==10
    plot( x, y_tp, ':r' ); 
  else
    plot( x, y_tp, '-r' );
  end
end

ps = linspace( 0, 1, 100 );
ys = 2 * exp( -60 * ps );
plot( ps, ys, '-g' );
xlabel('x'); ylabel('y(x)'); title('Ex 6.2.ii');

legend( 'shooting method', 'twopoint (n=10)', 'twopoint (n=50)', 'exact', 'location', 'best' );

%saveas( gcf, '../../WriteUp/Graphics/Chapter6/c_6_p_3_plot_ii.eps', 'epsc' );