function stiff_demo_2(step_range,Method)
%
% Implements the forward and backward Euler methods on 
% the stiff IVP 
%
%       y' = -1000*(y-t^2) + 2*t, y(0) = 0
%
% for a range of n-values and plots the log of the maximum 
% error vs. n.
%
% Inputs (optional): 
%  step_range = integer 2-vector [nlow, nhigh] specifying 
%                 the range of n-values. 
%  Method     = 'forward Euler', 'backward Euler', or 'both' (string). 
%
% Interesting exercise: Plot and compare the log maximum errors for 
% the two methods over various step ranges. (For example, use 
% step_range = [10 1000] or step_range = [480 520].) 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. 
%

% Set the defaults as needed and check the inputs. 

if nargin < 2, Method = 'both'; end
if nargin < 1, step_range = [480, 520]; end

if (step_range(1) > step_range(2))|(fix(step_range(1)) ~= step_range(1))|...
  (fix(step_range(2)) ~= step_range(2)), 
  error('The components of step_range must be integers in increasing order.');
end 

if ~(strcmp(Method,'forward Euler')|strcmp(Method,'backward Euler')|...
     strcmp(Method,'both')), 
  error('Method must be forward Euler, backward Euler, or both.'); 
end

% Complete the initialization. 

nlow = step_range(1); 
nhigh = step_range(2); 
nvals = [nlow:1:nhigh]';

errvals_fwd =zeros(size(nvals));
errvals_bwd =zeros(size(nvals));

% Loop over the range of steps. 

for nsteps = nlow:nhigh

  h = nsteps\1; 
  t = [0:h:1]';

% Implement the ODE solver(s). 

% Forward Euler. 

  if strcmp(Method,'forward Euler')|strcmp(Method,'both'),
    yvals = zeros(nsteps + 1,1); 
    for i = 1:nsteps 
      yvals(i+1,1) = yvals(i,1) + h*(-1000*(yvals(i,1)-t(i,1)^2) + 2*t(i,1));
    end
    errvals_fwd(nsteps-nlow+1,1) = log10(max(abs(yvals-t.^2)')); 
  end

% Backward Euler.

  if strcmp(Method,'backward Euler')|strcmp(Method,'both'),
    yvals = zeros(nsteps + 1,1); 
    for i = 1:nsteps 
      yvals(i+1,1) = (1+1000*h)\(yvals(i,1) + h*t(i+1,1)*(1000*t(i+1,1)+2)); 
    end
    errvals_bwd(nsteps-nlow+1,1) = log10(max(abs(yvals-t.^2)')); 
  end
end

% Plot the results. 

if strcmp(Method,'forward Euler'), 
  rlow = min(errvals_fwd)-1; 
  rhigh = max(errvals_fwd)+1;
  plot(nvals,errvals_fwd); axis([nlow nhigh rlow rhigh]);
  legend('forward Euler'); 
elseif strcmp(Method,'backward Euler'), 
  rlow = min(errvals_bwd)-1; 
  rhigh = max(errvals_bwd)+1;
  plot(nvals,errvals_bwd); axis([nlow nhigh rlow rhigh]);
  legend('backward Euler'); 
else 
  rlow = min([errvals_fwd;errvals_bwd])-1; 
  rhigh = max([errvals_fwd;errvals_bwd])+1;
  plot(nvals,errvals_fwd,nvals,errvals_bwd); axis([nlow nhigh rlow rhigh]);
  legend('forward Euler','backward Euler'); 
end
title('IVP dy/dt = -1000(y - t^2) + 2t, y(0)=0'); 
xlabel('Number of Steps'); ylabel('Log(Max |Error|)'); 
figure(1); drawnow; 
