%
%   NDEMO2.M   Solve y' = exp(-2x)-y with y(0) = 1  using Improved Euler Method
%
%   This is the Matlab version of the class demo used to illustrate the basics
%   of the Improved Euler method.  It calls the IMPEULER function file to 
%   actually implement the method.  The solution vector is plotted along with 
%   the exact solution to demonstrate the improved accuracy obtained with 
%   smaller step sizes (increased N).
%
%   Some additional computations are also performed to extract intermediate results 
%   for comparison to hand calculations.
%
%   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 = 0.5;  yo = 1; 
%
%   loop over different number of intervals (N)
      N = input('Input number of intervals, N, over the range [0 0.5] for EQN1: ');
%
      while N ~= 0
%   get numerical solution
        [x,y] = impeuler('eqn1',xo,xf,yo,N);
%   get function evaluations for comparison to hand calcs
        h = (xf-xo)/N;  k1 = h*feval('eqn1',x,y);  k2 = h*feval('eqn1',x+h,y+k1);
%   get exact solution for this problem (at same x values)
        ye = 2*exp(-x) - exp(-2*x);
%   compute max percent error 
        err = (y-ye)./ye;   errm = 100*max(abs(err));
%   echo numerical results to screen
        disp(' ')
        disp(['Improved EULER Results for N = ',num2str(N)])
        disp('   i     x-value  y-impeuler     k1        k2      y-exact     % error')
        for i = 1:N+1
        fprintf(1,'  %2i     %5.3f     %7.4f    %7.4f   %7.4f   %7.4f   %11.3e \n',...
               i,x(i),y(i),k1(i),k2(i),ye(i),err(i)*100)
        end

        disp(' ')
%   plot solution
        nfig = nfig+1;   figure(nfig)
        plot(x,y,'r.',x,ye,'g-','LineWidth',2),grid
        title(['Comparison of Solutions to y'' = exp(-2x) - y (N = ',num2str(N),')'])
        xlabel('x value'),ylabel('y-value')
        legend('ImpEuler Method','Exact Soln')
        r = axis;   r(2) = xf;   axis(r);
        xt = r(1) + 0.05*(r(2)-r(1));   yt = r(3) + 0.05*(r(4)-r(3));
        text(xt,yt,['Max Error = ',num2str(errm),' %'])
%   try another N?
        N = input('Try a different N (enter zero to quit)? ');
      end
%
%   end of demo
%
