%  
%   HT1D_SS2.M   Heat Transfer in a 1-D Finite Bar Using the STATE SPACE Method
%                   (1-D Space-Time Heat Transfer Example #2 in Class Notes)
%  
%   State Space 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 
%
%   Uses function file HT1D_SS2A.M within call to ODE23   
%
%   File prepared by J. R. White, UMass-Lowell (updated Dec. 2007)
%

%
%   getting started
      clear all,  close all,  nfig = 0;
      global A b 
%
%   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 increments and time domain parameters
      Nx = 24;    dx = L/(Nx+1);     x = dx:dx:L-dx;    % spatial grid
      to = 0;     tf = 40*60;                           % time interval (sec)
%
%   set edit parameters  and initial temperature distribution     
      Ne = 6;                 % number of edit intervals (including initial time)
      DT = (tf-to)/Ne;        % interval width
      t = linspace(to,tf,Ne); % times of interest for spatial snapshots
      f = u0*ones(size(x'));  % initial temperature distribution
%
%   setup coeff matrix
      A = zeros(Nx,Nx);   b = zeros(Nx,1);    aa = alf/(dx*dx);
      b(1) = aa*ul;   A(1,1) = -2*aa;  A(1,2) = aa;       % left boundary 
      for i = 2:Nx-1
        A(i,i-1) = aa;   A(i,i) = -2*aa;   A(i,i+1) = aa; % interior nodes 
      end
      b(Nx) = aa*ur;  A(Nx,Nx-1) = aa; A(Nx,Nx) = -2*aa;  % right boundary
%
%   now use ode23 to solve system of ODEs -- gives spatial profile at each time point
      [tt,u] = ode23('ht1d_ss2a',t,f);
      u = u';  % store solution as u(x,t1), u(x,t2), ... where i^th column is for t_i
%
%   add boundary points to solution for ploting
      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\_SS2:  Bar Temperature Profiles at Several Times')                  
      grid,xlabel('Distance (m)'),ylabel('Temperature (C)')       
      legend(char(st),'Location','NorthWest')
%
%   end of problem
                                                                 

