%                                                                             
%   FD.M   Example of Finite Difference Method for Solving ODEs          
%                                                                             
%   This demo uses the finite difference technique as the solution method                 
%   for a particular ODE.  The mesh interval is constant.  Two different            
%   approximations for the first derivative are used.  This comparison 
%   demonstrates that the best approximation for the derivative terms 
%   leads to a more efficient numerical scheme.  One should always use 
%   the best approximation available for the given problem!  
%
%   The sensitivity to mesh size is also addressed by looping over four
%   different values of N (number of unknowns in the problem).        
%                                                                             
%   Description of given differential equation                                
%     y'' + 3xy' + 7y  =  cos(2x)      with  y(0) = 1  and  y(pi) = 0         
% 
%   A related Matlab file, FDS.M, is also available (this contains a simple 
%   version of this example - Case A only - with no mesh sensitivity study).          
%  
%   File prepared by J. R. White, UMass-Lowell (Aug. 2003).
%

%
%   getting started
      clear all, close all, nfig = 0;
%                                                                           
%   set domain limits and boundary conditions
      xo = 0;   xf = pi;  yxo = 1;   yxf = 0; 
%
%   loop over four values of N to evaluate sensitivity to mesh size
      NN = [20 40 80 160];   
%
      for n = 1:4
        N = NN(n);
%
%   compute interval size and discrete x vector
        dx = (xf-xo)/(N+1);   dx2 = dx*dx;    x = (xo+dx):dx:(xf-dx);
%                                                                             
%                                                                             
%   Approximations:  y'' = [y(i-1) - 2y(i) + y(i+1)]/(dx^2)  error => (dx^2)
%     (Case A)       y'  = [y(i+1) - y(i-1)]/(2dx)           error => (dx^2)     
%                                                                             
%                                                                             
%   setup matrix eqns. (treat boundary terms as special cases)  
        a = zeros(N,N);   b = zeros(N,1);
        for i = 2:N-1                                                          
          a(i,i-1) = 1-3*x(i)*dx/2;    
          a(i,i) = -2+7*dx2;                                                     
          a(i,i+1) = 1+3*x(i)*dx/2;     
          b(i) = dx2*cos(2*x(i));     
        end                                                                       
%   left boundary                                                             
        a(1,1) = -2+7*dx2;                                                        
        a(1,2) = 1+3*x(1)*dx/2;  
        b(1) = dx2*cos(2*x(1))-(1-3*x(1)*dx/2)*yxo;
%   right boundary                                                            
        a(N,N-1) = 1-3*x(N)*dx/2;                                          
        a(N,N) = -2+7*dx2;                                                  
        b(N) = dx2*cos(2*x(N))-(1+3*x(N)*dx/2)*yxf;                         
%   solve system of eqns                                                      
        y = a\b;                                                                 
%   add boundary points to solution for plotting                              
        za = [yxo y' yxf];   xa = [xo x xf];
%                                                                               
%                                                                             
%   Approximations: y'' = [y(i-1) - 2y(i) + y(i+1)]/(dx^2)  error => (dx^2) 
%     (Case B)      y'  = [y(i+1) - y(i)]/(dx)              error => (dx)           
%                                                                             
%                                                                             
%   setup matrix eqns. (treat boundary terms as special cases)  
        a = zeros(N,N);   b = zeros(N,1);
        for i = 2:N-1                                                          
          a(i,i-1) = 1;    
          a(i,i) = -(2+3*x(i)*dx-7*dx2); 
          a(i,i+1) = 1+3*x(i)*dx;     
          b(i) = dx2*cos(2*x(i));     
        end                                                                       
%   left boundary                                                             
        a(1,1) = -(2+3*x(1)*dx-7*dx2);
        a(1,2) = 1+3*x(1)*dx;  
        b(1) = dx2*cos(2*x(1))-yxo;
%   right boundary                                                            
        a(N,N-1) = 1;                                          
        a(N,N) = -(2+3*x(N)*dx-7*dx2);
        b(N) = dx2*cos(2*x(N))-(1+3*x(N)*dx)*yxf;                         
%   solve system of eqns                                                      
        y = a\b;                                                                 
%   add boundary points to solution for plotting                              
        zb = [yxo y' yxf];   xb = [xo x xf];
%                                                                             
%   plot results for both cases                                                             
        t = '220+n';  
        subplot(eval(t)),plot(xa,za,'r-',xb,zb,'g--','LineWidth',2) 
        axis([0 3.2 -2 1]);
        title(['FD:  FD Method (',num2str(N),' pts)']) 
        xlabel('x values'),ylabel('y values'),grid                                
        legend('Case A','Case B')                                           
%
      end
%
%   end of demo
