%
%   GFSERIES1.M   Demo using Eigenfunctions of a Sturm-Liouville
%               for a Generalized Fourier Series Representation for f(x) 
%
%   This problem follows an example given in class.  There are a number of 
%   steps involved, as follows:
%     1. use fzero to find the roots to the eigencondition equation (this 
%        determines the eigenvalues for the problem).
%     2. use quadl to numericallly integrate the orthogonality statement
%        to prove that the eigenfucntions are indeed orthogonal and to 
%        determine the norm of the functions (when m = n)
%     3. actually expand two different functions, f1(x) and f2(x), using 
%        the orthogonal eigenfunctions (these will also need quadl to find 
%        the expansion coeffs via numerical integration).
%     4. plot the eigenfunction expansions along with the exact functions,
%        f1(x) and f2(x), to see how well the truncated generalized Fourier 
%        series matches the exact expressions.
%
%   This should be instructive -- let's give it a try...
%
%   File prepared by J. R. White, UMass-Lowell (Nov. 2003)
%

%
%   getting started
      clear all, close all, nfig = 0;
      global h L ALPM ALPN
%
%   set problem parameters
      L = 2.5;  h = 2;      
      Nx = 251;   x = linspace(0,L,Nx)';
%
%   Part 1  --  Find eigenvalues
      Ne = 25;      % number of eigenvlaues (also max order of expansion)
      alp = zeros(Ne,1);
      for n = 1:Ne
        g1 = (2*n-1)*pi/2;  g2 = (2*n+1)*pi/2;
        pn = fzero('gfseries1a',[g1 g2]);
        alp(n) = pn/L;
      end
%
%   Part 2  --  Evaluate the orthogonality condition
      ymyn = zeros(Ne,Ne);
      for n = 1:Ne
        ALPN = alp(n);
        for m = 1:Ne
          ALPM = alp(m);
          ymyn(m,n) = quadl('gfseries1b',0,L);
        end
      end
%
%     ------------------------------------------------
%     |       Case 1  --  f1(x) = 3*exp(-2x)         |
%     ------------------------------------------------
%
%   Part 3a  --  Find coeff for f1(x) = 3*exp(-2x)
      a1 = zeros(Ne,1);
      for n = 1:Ne
        ALPN = alp(n);
        a1(n) = quadl('gfseries1c1',0,L)/ymyn(n,n);
      end
%
%   Part 3b  --  Calc exact function 
      f1e = 3*exp(-2*x);
%
%   Part 3c  --  Calc Fourier Series approx
      f1a = zeros(Nx,1);   tol = 0.01;   mrerr = 1.0;   n = 0;
      while mrerr > tol   &    n < Ne
        n = n+1;     
        ff = a1(n)*sin(alp(n)*x);   f1a = f1a + ff;
        rerr = ff(2:Nx-1)./f1a(2:Nx-1);   mrerr = max(abs(rerr));
      end
%
%   Part 4  --  Plot curves
      nfig = nfig+1;   figure(nfig)
      plot(x,f1e,'g-',x,f1a,'r--','LineWidth',2),grid
      title('GFSeries1:  Accuracy of Truncated Fourier Series for f1(x) = 3*exp(-2*x)')
      ylabel('f1(x)'), xlabel('x value')
      legend('Exact f1(x)','Approx f1(x)')
      if n == Ne,  gtext(['not converged with ',num2str(Ne),' terms']);  end
      if n < Ne,   gtext([num2str(n),' terms for convergence']);  end
%
%     ---------------------------------------------------------
%     |       Case 2  --  f2(x) = 3*exp(-2*x) + x - 3         |
%     ---------------------------------------------------------
%
%   Part 3a  --  Find coeff for f2(x) = 3*exp(-2*x) + x - 3
      a2 = zeros(Ne,1);
      for n = 1:Ne
        ALPN = alp(n);
        a2(n) = quadl('gfseries1c2',0,L)/ymyn(n,n);
      end
%
%   Part 3b  --  Calc exact function 
      f2e = 3*exp(-2*x) + x - 3;
%
%   Part 3c  --  Calc Fourier Series approx
      f2a = zeros(Nx,1);   tol = 0.01;   mrerr = 1.0;   n = 0;
      while mrerr > tol   &    n < Ne
        n = n+1;     
        ff = a2(n)*sin(alp(n)*x);   f2a = f2a + ff;
        rerr = ff(2:Nx-1)./f2a(2:Nx-1);   mrerr = max(abs(rerr));
      end
%
%   Part 4  --  Plot curves
      nfig = nfig+1;   figure(nfig)
      plot(x,f2e,'g-',x,f2a,'r--','LineWidth',2),grid
      title('GFSeries1:  Accuracy of Truncated Fourier Series for f2(x) = 3*exp(-2*x) + x - 3')
      ylabel('f2(x)'), xlabel('x value')
      legend('Exact f2(x)','Approx f2(x)')
      if n == Ne,  gtext(['not converged with ',num2str(Ne),' terms']);  end
      if n < Ne,   gtext([num2str(n),' terms for convergence']);  end
%
%   end of demo