%  
%   HT1D_EX2.M   Heat Transfer in a 1-D Finite Bar Using EXPLICIT Numerical Method
%                   (1-D Space-Time Heat Transfer Example #2 in Class Notes)
%  
%   Explict Numerical Solution to the following problem: 
%        ut(x,t) = alf*uxx(x,t)  
%   with u(0,t) = ul   u(L,t) = ur   u(x,0) = u0 
%
%   File generated by J.R. White, UMass-Lowell (updated Dec. 2007)
%
   
%
%   getting started
      clear all,  close all,  nfig = 0;
%
%   problem data                                                         
      alf = .0001;     % m^2/s  thermal diffusivity                            
      ul = 100;        % C      fixed temp at left end                         
      ur = 400;        % C      fixed temp at right end                        
      u0 = 70;         % C      initial uniform temp of bar                    
      L = 1;           % m      length of bar 
%
%   set dx and dt increments
      Nx = 24;    dx = L/(Nx+1);     x = dx:dx:L-dx; 
      Fo = input('Input value of Fourier Modulus (Fo): ');
      dt = Fo*dx*dx/alf;
      tf = 40;         % approximate final time for simulation (min.)
      tf = tf*60;   Nt = ceil(tf/dt);    
      fprintf(1,'\nDiscrete Values for this run: \n');
      fprintf(1,'Time interval size (sec):      %6.3f \n',dt);
      fprintf(1,'Spatial interval size (m):     %6.3f \n',dx);
      fprintf(1,'Number of time intervals (Nt): %6i   \n',Nt);  
      Ne = 5;   Nte = ceil(Nt/Ne);
      fprintf(1,'Number of edit intervals (Ne): %6i   \n',Ne);  
      fprintf(1,'Number of time steps per edit interval (Nte): %6i \n',Nte);
      Ne = Ne+1;       % number of edit points (add one for initial state)
%
%   setup storage for temperature profiles (one column for each edit time point)                                                      
      u = zeros(Nx,Ne);   t = zeros(1,Ne);
      u(:,1) = u0*ones(size(x))';   t(1) = 0;   uold = u(:,1);  
%
%   setup coeff matrix (ok for small problems, not very efficient in general)
      A = zeros(Nx,Nx);   b = zeros(Nx,1);
      b(1) = Fo*ul;    A(1,1) = 1-2*Fo;   A(1,2) = Fo;       % left boundary 
      for i = 2:Nx-1
        A(i,i-1) = Fo;   A(i,i) = 1-2*Fo;   A(i,i+1) = Fo;   % interior nodes 
      end
      b(Nx) = Fo*ur;   A(Nx,Nx-1) = Fo;   A(Nx,Nx) = 1-2*Fo; % right boundary
%
%   do loop over time
      tt = 0;
      for j = 2:Ne
        nc = 0;
        while nc < Nte
          nc = nc+1;   unew = A*uold + b;   uold = unew;   tt = tt+dt;
        end
        t(j) = tt;    u(:,j) = unew;
      end
%
%   add boundary points to solution for plotting
      xp = [0 x L];   up = zeros(Nx+2,Ne);
      up(:,1) = [u0; u(:,1); u0];
      for j = 2:Ne;    up(:,j) = [ul; u(:,j); ur];   end
%
%   set color and marker code for creating plots
      Ncm = 6;  scm = ['r-';'g:';'b-';'m:';'c-';'k:']; 
%
%   plot curves of up for various times
      st = zeros(Ne,9);
      nfig = nfig+1;   figure(nfig) 
      jj = 0;                       
      for j = 1:Ne
        jj = jj+1;  if jj > 6, jj = 1; end
        plot(xp,up(:,j),scm(jj,:),'LineWidth',2); hold on 
        st(j,:) = sprintf('%5.1f min',t(j)/60);
      end
      axis([0 L 0 400]); 
      hold off                          
      title(['HT1D\_EX2:  Bar Temperature Profile at Several Times (Fo = ',...
             num2str(Fo),')'])                  
      grid,xlabel('Distance (m)'),ylabel('Temperature (C)')       
      legend(char(st),'Location','NorthWest')
%
%   end of problem
                                                                 

