function stiff_demo_1(nsteps,Method)
%
% Implements the forward and backward Euler methods on 
% the stiff IVP 
%
%       y' = -1000*(y-t^2) + 2*t, y(0) = 0
%
% Inputs: 
%  nsteps = number of steps across [0,1]; may be a scalar or a vector. 
%  Method = 'forward Euler' or 'backward Euler' (string). 
%
% Interesting exercise: Compare the maximum errors for the two methods 
% with various values of nsteps. In particular, note 
%  * the huge errors with forward Euler with 0 < nsteps < 500, 
%  * the dramatic onset of stability of forward Euler around nstep = 500, 
%  * the steady O(h) behavior of the backward Euler errors for all 
%    values of h = 1/nsteps. 
%

% Loop over the number of steps. 

for k = 1:length(nsteps), 

% Initialize.
 
  h = nsteps(k)\1; 
  t = [0:h:1]';
  yvals = zeros(nsteps(k) + 1,1); 

% Call the ODE solver and print the results. 

  if strcmp(Method,'forward Euler'), 
    for i = 1:nsteps(k) 
      yvals(i+1,1) = yvals(i,1) + h*(-1000*(yvals(i,1)-t(i,1)^2) + 2*t(i,1));
    end
  elseif strcmp(Method,'backward Euler'), 
    for i = 1:nsteps(k) 
      yvals(i+1,1) = (1+1000*h)\(yvals(i,1) + ... 
  		   h*t(i+1,1)*(1000*t(i+1,1)+2)); 
    end
  else 
    error('Method must be either forward Euler or backward Euler.')
  end
  maxerr = max(abs(yvals-t.^2)'); 
  fprintf('\n %s: nsteps = %d    max. error %e ', Method, nsteps(k), maxerr); 
end

fprintf('\n\n');
