%
%   PSERIES_EX1.M   Evaluate the solution to a Linear ODE using
%                     Power Series Expansion Method
%
%   This demo evaluates the infinite series expansion that was
%   obtained using the Power Series Expansion technique when solving
%   the following IVP:
%      (x^2 - 4)y'' + 3xy' + y = 0
%   with  y(0) = 4  and  y'(0) = 1
%
%   It illustrates how to evaluate the recurrence relation and how 
%   to check for convergence of the series. 
%
%   As a check, it also uses ode23 to solve the given IVP --  which 
%   allows comparison of the analytical (Power Series) and numerical 
%   solutions.
%
%   File prepared by J. R. White, UMass-Lowell (Oct. 2003)
%

%
%   getting started
      clear all,  close all,   nfig = 0;
%
%   define range of x values and initial conditions
      xo = 0;  xf = 1.8;  yo = 4;  ypo = 1;        
%
%   evaluate numerical solution with Matlab's built-in ode solver
      zo = [yo ypo]';
      [xn,z] = ode23('pseries_ex1a',[xo xf],zo);
      yn = z(:,1);                          
%
%   define x values and max number of terms for analytical solution
      Nx = 30;  x = linspace(xo,xf,Nx)';
      maxT = 100;  tol = 0.0001;     
%
%   evaluate y1(x) solution using infinite series expansion
      n = 0;  mrerr = 1.0;
      a0 = 4;   T = a0*ones(size(x));  y1 = T; 
      while mrerr > tol  &  n < maxT
        n = n+1;   m = 2*n-2; 
        r = x.*x*(m+1)/(4*(m+2));   T = r.*T;
        y1 = y1 + T;
        ii = find(y1);  rerr = T(ii)./y1(ii);   mrerr = max(abs(rerr));
      end
      NT1 = n;
%
%   evaluate y2(x) solution using infinite series expansion 
      n = 0;  mrerr = 1.0;
      a1 = 1;   T = a1*x;  y2 = T; 
      while mrerr > tol  &  n < maxT 
        n = n+1;   m = 2*n-1;
        r = x.*x*(m+1)/(4*(m+2));   T = r.*T;
        y2 = y2 + T;
        ii = find(y2);  rerr = T(ii)./y2(ii);   mrerr = max(abs(rerr));
      end
      NT2 = n;
%
%   put analytical solution together
      y = y1 + y2;
%
%   plot functions
      nfig = nfig+1;   figure(nfig)
      plot(x,y,'r--',xn,yn,'bo','LineWidth',2), grid
      title('PSeries\_Ex1:  Power Series Solution  --  Example #1'); 
      xlabel('x-values'),ylabel('y-values')
      legend(['Power Series (NT1 = ',num2str(NT1),' & NT2 = ',num2str(NT2),')'], ...
             'ODE23 Numerical Solution');
%
%   end of program