%
%                 INODE - Interactive Numerical solution of ODEs
%
%   This M-file solves the chosen differential equation using one of four  
%   numerical integration schemes -- Euler, Improved Euler, 4th Order RK, or the 
%   adaptive predictor-corrector method within Matlab - ODE23.
%
%   The function file containing the definition of f(x,y) is selected interactively
%   (must be available a priori) or the user can actually type in the function
%   f(x,y).  The second option calls a generic function file, sfode.m, which 
%   evaulates the function f(x,y) input via keyboard.  The integration limit 
%   are also set by the user. 
%
%   This file makes it pretty easy to solve a large variety of first order ODEs.
%
%   File prepared by J. R. White, UMass-Lowell (original March 1998)
%     --> modified slightly (Oct. 1998)
%

%
%   getting started
      if exist('nfig') == 0,  nfig = 0;  end   % initialize nfig if it does not exist
      global sfunc
%
%   identify equation to be solved
      fopt = menu('Select option for identifying ODE to solve:', ...
                  'Function file containing ODE definition already exists', ...
                  'Enter function f(x,y) from keyboard interactively');
      if fopt == 1
        [name]=uigetfile('eqn*.m'); name = lower(name);
        nn = findstr(name,'.m');    fxy = name(1:nn-1); 
        ftitl = fxy;
      else
        sfunc = input('  Enter desired f(x,y): ','s');
        ftitl = sfunc;
        fxy = 'sfode';
      end
%
%   define some initial parameters
      disp(['  Identify range of independent variable and initial conditions for ', ...
           ftitl,': ']);
      xo = input('    Input initial x-value: ');
      yo = input('    Input initial y-value: ');
      xf = input('    Input final x-value: ');
      disp(' ')
%
%   select technique to use
      meth = menu('Select methd to use', ...
                  ' EULER method (N = 100)', ...
                  ' IMPEULER method (N = 100)', ...
                  ' RK4 method (N = 100)', ...
                  ' ODE23 adaptive step control (tol = 1e-3)');
%   now get actual solutions
      N = 100;   tol = 1e-3;   options = odeset('RelTol',tol);
      if meth == 1,  [x,y] = euler(fxy,xo,xf,yo,N);          end
      if meth == 2,  [x,y] = impeuler(fxy,xo,xf,yo,N);       end
      if meth == 3,  [x,y] = rk4(fxy,xo,xf,yo,N);            end
      if meth == 4,  [x,y] = ode23(fxy,[xo xf],yo,options);  end
%        
%   plot solutions
      nfig = nfig+1;   figure(nfig)
      plot(x,y),grid
      title(['Numerical Solution Techniques for y'' = ', ...
             ftitl,' (Method #',num2str(meth),')'])
      xlabel('x value'),ylabel('y-value')
      r = axis;   r(1) = xo;  r(2) = xf;   axis(r);
%
%   end of demo
%
