%                                                                             
%   FDS.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.  
%
%   Description of given differential equation                                
%     y'' + 3xy' + 7y  =  cos(2x)      with  y(0) = 1  and  y(pi) = 0         
% 
%   A related Matlab file, FD.M, is also available (this contains a more detailed 
%   version of this example - compares two different approximations for y' and
%   performs a 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; 
%
%   user defined number of unknowns in problem 
      N = input('Input number of unknowns in problem (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];
%                                                                             
%   plot results 
      nfig = nfig+1;  figure(nfig)                                                             
      plot(xa,za,'r-','LineWidth',2) 
      axis([0 3.2 -2 1]);
      title(['FDs:  FD Method for Case A (',num2str(N),' pts)']) 
      xlabel('x values'),ylabel('y values'),grid                                
%
%   end of demo
