% Sample driver for Problem 2
% You need to implement nonlinshoot.m and update problem data

% 11.2.3 a
fprintf('---------------- 11.2.3 a ----------------\n');
f   = @(t,y,yp) y^3-y*yp;
fy  = @(t,y,yp) 3*y^2 - yp;
fyp = @(t,y,yp) -y;
ytrue = @(t) 1/(t+1);
a=1; b=2; alpha=1/2; beta=1/3; n=5; h=(b-a)/n;
maxit=100; tol=1e-4;
y = nonlinshoot(a,b,n,alpha,beta,f,fy,fyp,tol,maxit);
fprintf('%3s %11s %11s\n','t','y','y true');
for i=1:n+1,
 ti=a+(i-1)*h;
 fprintf('%3.1f %11.8f %11.8f\n',ti,y(i),ytrue(ti));
end;

% 11.2.3 b
fprintf('---------------- 11.2.3 b ----------------\n');
f   = @(t,y,yp) 2*y^3 - 6*y -2*t^3;
fy  = @(t,y,yp) 6*y^2 -6;
fyp = @(t,y,yp) 0;
ytrue = @(t) t + 1/t;
maxit=100; tol=1e-4;
a=1; b=2; alpha=2; beta=5/2; n=5; h=(b-a)/n;
y = nonlinshoot(a,b,n,alpha,beta,f,fy,fyp,tol,maxit);
fprintf('%3s %11s %11s\n','t','y','y true');
for i=1:n+1,
 ti=a+(i-1)*h;
 fprintf('%3.1f %11.8f %11.8f\n',ti,y(i),ytrue(ti));
end;

% typical output
%
%  >> prob2_sample
%  ---------------- 11.2.3 a ----------------
%  it=  0, z=-0.1666667, phi=  7.4898e-02
%  it=  1, z=-0.2503546, phi= -3.1549e-04
%  it=  2, z=-0.2500051, phi= -6.3503e-09
%    t           y      y true
%  1.0  0.50000000  0.50000000
%  1.2  0.45454618  0.45454545
%  1.4  0.41666752  0.41666667
%  1.6  0.38461607  0.38461538
%  1.8  0.35714323  0.35714286
%  2.0  0.33333333  0.33333333
%  ---------------- 11.2.3 b ----------------
%  it=  0, z= 0.5000000, phi=  4.0955e+01
%  it=  1, z= 0.3842470, phi=  1.6131e+01
%  it=  2, z= 0.2523326, phi=  6.0554e+00
%  it=  3, z= 0.1198815, phi=  1.9112e+00
%  it=  4, z= 0.0289893, phi=  3.6810e-01
%  it=  5, z= 0.0020521, phi=  2.0387e-02
%  it=  6, z= 0.0003799, phi=  7.0127e-05
%    t           y      y true
%  1.0  2.00000000  2.00000000
%  1.2  2.03369066  2.03333333
%  1.4  2.11466252  2.11428571
%  1.6  2.22530754  2.22500000
%  1.8  2.35576849  2.35555556
%  2.0  2.50007013  2.50000000

