%
%   ODENEX3.M   Solution of ODEs with MATLAB  (IVP)  
%
%   This demo gives an example of how to use the ODE23 routine to integrate
%   a simple first order differential equation. A slight extension to vectors
%   allows the solution of systems of first order ODEs.  This extension will be 
%   treated in later examples (see section on General IVPs).
%
%   This example also illustrates how to evaluate and plot a function in implicit
%   form.  This requires use of Matlab's FZERO root finding routine.
%
%   also see ODENEX3A.M - function file for use in FZERO
%            ODENEX3B.M - function file for use in ODE23
%
%   This example solves   y' = 1/(x + y*exp(y))   with    y(0) = 1
%   and compares to exact soln   -x*exp(-y) + 0.5*y^2 = 0.5   (implicit form)
%
%   File prepared by J. R. White, UMass-Lowell (Sept 2006)
%

%
%   getting started
      clear all, close all, nfig = 0;
%
%   initial setup 
      xo = 0;   xf = 3;   yo = 1;   tol = .00001;
      global XX
%
%   evaluate exact solution
%   The implicit form can be written as f(x,y) = -x*exp(-y) + 0.5*y^2 - 0.5 = 0 
%   Thus, to plot y versus x, we find the 'roots' of this implicit eqn 
%   for each value of x, and then plot the results.
%
%   first set the range for x
      Nx = 31;   xe = linspace(xo,xf,Nx);   ye = zeros(size(xe));   
%
%   now evaluate the implicit eqn at each point (uses previous y-value as guess)
      xe(1)= xo;    ye(1) = yo;
      for i = 2:Nx
        XX = xe(i);
        ye(i) = fzero('odenex3a',ye(i-1));  
      end
%   
%   now evaluate the numerical solution (using ODE23)
      options = odeset('RelTol',tol);
      [x1,y1] = ode23('odenex3b',[xo xf],yo,options);
%
%   now plot results (numerical(points), exact(solid))
      nfig = nfig+1;  figure(nfig)
      plot(xe,ye,'b-',x1,y1,'ro','LineWidth',2)
      title('ODEnex3:  Numerical & Analytical Solutions to 1st Order ODE')
      xlabel('x values'),ylabel('y(x)'),grid
      legend('EXACT','ODE23');
%
%   end of simulation