function [] = prob_2_10_pt2
% EX_2_10_PT2 -
%
% Solve y' = f(t,y) with the trapezoidal rule:
%
% ynp1 = y_n + 0.5 * h * ( f(tn,yn) + f(tnp1,ynp1) );
% 
% using functional iterations.
% 
% 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);
  ynp1 = AM2si(tn,yn,h,@ode);
  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 [ynp1] = AM2si(tn,yn,h,f)
% AM2SI - 
% 
% 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.
% 
%-----

nIters = 1;
ynp1 = yn + 0.5 * h * ( feval(f,tn,yn) + feval(f,tn+h,yn) );
while( (norm(ynp1-yn) >= 1e-3 * norm(ynp1)) && nIters < 11 )
  nIters = nIters+1;
  ynp1 = yn + 0.5 * h * ( feval(f,tn,yn) + feval(f,tn+h,ynp1) );
end

return;