function [] = prob_2_10_pt1
% EX_2_10 - 
% 
% 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 = Heun(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] = Heun(tn,yn,h,f)
% HEUN - 
% 
% 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.
% 
%-----

ync1 = yn + h * feval(f,tn,yn);;
ynp1 = yn + 0.5 * h * ( feval(f,tn,yn) + feval(f,tn+h,ync1) );
return;