%
%   NDEMO4.M   Solve y' = exp(-2x)-y with y(0) = 1  using a variety of method
%
%   This M-file solves the class demo using a full selection of numerical 
%   integration schemes -- Euler, Improved Euler, 4th Order RK, and the adaptive
%   predictor-corrector method within Matlab - ODE23.
%   
%   The solution vectors are plotted along with the exact solution to demonstrate 
%   the accuracy obtained with the different methods.
%
%   The function file EQN1.M containing the function definition is also required.
%
%   File prepared by J. R. White, UMass-Lowell (March 2006)
%

%
%   getting started
      clear all;   close all;   nfig = 0;
%
%   define some initial parameters
      xo = 0;   xf = 5.0;  yo = 1; 
%
%   get exact solution for this problem 
      xe = linspace(xo,xf,51);   ye = 2*exp(-xe) - exp(-2*xe);
%
%   get numerical solution using adpative step control technique (ODE23)
      tol = 1e-4;   options = odeset('RelTol',tol);
      [x4,y4] = ode23('eqn1',[xo xf],yo,options);
%
%   get numerical solutions for fixed step size
      N = input('Input number of intervals, N, over the range [0 5.0] for EQN1: ');
      while N ~= 0
        [x1,y1] = euler('eqn1',xo,xf,yo,N);
        [x2,y2] = impeuler('eqn1',xo,xf,yo,N);
        [x3,y3] = rk4('eqn1',xo,xf,yo,N);
%   plot solution
        nfig = nfig+1;   figure(nfig)
        plot(xe,ye,'r-',x1,y1,'g.',x2,y2,'b+',x3,y3,'ms',x4,y4,'co','LineWidth',2),grid
        title(['Comparison of Solutions to y'' = exp(-2x) - y (N = ',num2str(N),')'])
        xlabel('x value'),ylabel('y-value')
        legend('Exact Soln','Euler Method','ImpEuler Method','RK4 Method','ODE23 Method')
        r = axis;   r(2) = xf;   axis(r);
%   try another N?
        N = input('Try a different N (enter zero to quit)? ');
      end
%
%   end of demo
%
