%function prob_7_2()
% 
% Written by:
% -- 
% John L. Weatherwax                2006-05-21
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

r = 0.06; 
tFinal = 12; 
h=1/12; 

y0=100; 

% Euler's method compounding yearly ...
%
h=1;
nSteps=floor(tFinal/h);
ys1=zeros(1,nSteps+1);
ys1(1)=y0;
for si=1:nSteps,
  ynp1 = ys1(si) + h*r*ys1(si);
  ys1(si+1) = ynp1; 
end

% Euler's method compounding monthly ...
%
h=1/12;
nSteps=floor(tFinal/h); 
ys2=zeros(1,nSteps+1);
ys2(1)=y0; 
for si=1:nSteps,
  ynp1 = ys2(si) + h*r*ys2(si);
  ys2(si+1) = ynp1;
end

% Midpoint method compounding monthly ...
%
h=1/12;
nSteps=floor(tFinal/h); 
ys3=zeros(1,nSteps+1); 
ys3(1)=y0;
for si=1:nSteps,
  s1 = r*ys3(si); 
  s2 = r*(ys3(si)+0.5*h*s1);
  ys3(si+1) = ys3(si) + h*s2;
end

% Trapezoid method compounding monthly ...
%
h=1/12;
nSteps=floor(tFinal/h); 
ys4=zeros(1,nSteps+1); 
ys4(1)=y0;
for si=1:nSteps,
  s1 = r*ys4(si); 
  s2 = r*(ys4(si)+h*s1);
  ys4(si+1) = ys4(si) + 0.5*h*(s1+s2);
end

% BS23 algorithm compounding monthly ...
%
h=1/12;
nSteps=floor(tFinal/h); 
ys5=zeros(1,nSteps+1); 
ys5(1)=y0;
for si=1:nSteps,
  s1 = r*ys5(si); 
  s2 = r*(ys5(si)+0.50*h*s1);
  s3 = r*(ys5(si)+0.75*h*s2);
  ys5(si+1) = ys5(si) + (1/9)*h*(2*s1+3*s2+4*s3);
end

% Continious compounding ...
%
h=1/12;
ts=0:h:12; 
ys6=y0*exp(r*ts);


% Plot all results: 
% 
figure; plot( 0:1:12, ys1, '-x' ); grid on; hold on; 
Ys = [ ys2; ys3; ys4; ys5; ys6 ]; 
plot( ts, Ys );