function [] = prob_2_10_pt3
% EX_2_10_PT3 -
%
% Solve y' = f(t,y) with the trapezoidal rule:
%
% ynp1 = y_n + 0.5 * h * ( f(tn,yn) + f(tnp1,ynp1) );
% 
% using quasi-newtons method.
% 
% Written by:
% -- 
% John L. Weatherwax                2005-03-10
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

a = 0; b = 1; h = 0.1; 
nSteps = round((b-a)/h);

ys = [ 0; 1];
for ii=1:nSteps,
  tn = a + (ii-1)*h;
  yn = ys(:,end);
  [tnp1,ynp1] = AM2ni(tn,yn,h,@ode,@dode);
  ys = [ ys, ynp1 ];
end

ts = a + h*(0:nSteps);
yTruth = [ sin(ts); cos(ts) ];

figure; 
subplot(1,2,1);
plot( ts, ys(1,:), '-r', ts, yTruth(1,:), '-g' );
title( 'y1 approx & truth' );
subplot(1,2,2);
plot( ts, ys(2,:), '-r', ts, yTruth(2,:), '-g' );
title( 'y2 approx & truth' );

return;


%------------------------------------------------------------------------
%------------------------------------------------------------------------

function dydt = ode(t,y)
dydt = [ 0 1; -1 0 ] * y;

%------------------------------------------------------------------------
%------------------------------------------------------------------------

function J = dode(t,y)
J = [ 0 1; -1 0 ];


%------------------------------------------------------------------------
%------------------------------------------------------------------------

function [tnp1,ynp1] = AM2ni(tn,yn,h,f,dfdy)
% AM2NI - 
% 
% Written by:
% -- 
% John L. Weatherwax                2005-03-10
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

% Since this is the quasi-newton method we evaluate the
% derivative at yn only once and use it for all iterations:
J = feval(dfdy,tn,yn);
M = eye(size(J)) - 0.5*h*J;
[L,U]=lu(M);

% Initialize:
tnp1 = tn+h;
nIters = 0; ynp1 = yn;

% Do one iteration outside the while loop:
b = yn + 0.5 * h * feval(f,tn,yn) + 0.5 * h * feval(f,tnp1,ynp1) - ynp1;
deltam = U \ ( L \ b ); 
ynp1 = ynp1 + deltam;

% Do remaining iterations (if needed):
nIters = 1;
while( (norm(deltam) >= 1e-3 * norm(ynp1)) && nIters < 101 )
  b = yn + 0.5 * h * feval(f,tn,yn) + 0.5 * h * feval(f,tnp1,ynp1) - ynp1;
  deltam = U \ ( L \ b );
  ynp1 = ynp1 + deltam;
  nIters = nIters+1;
end

nIters

return;