%
%   HT1D_SOV1.M   Heat Transfer in a 1-D Finite Bar (Example 10.1 in Class Notes)
%
%   Analytical Solutionn using Separation of Variables (SOV) method for the 
%   following problem:
%     ut(x,t) = alf*uxx(x,t)   with u(0,t) = 0   u(L,t) = 0   u(x,0) = f(x)
%     where    u(x,0) = f(x) = x    if 0 < x <L/2       
%                             L-x  if L/2 < x < L 
%
%   The goal here is to evaluate and plot the temperature profile in the 
%   1-D finite bar of length L at various time points.  This is an easy
%   way to visualize the full space-time solution, u(x,t).
%
%   File prepared by J. R. White, UMass-Lowell (Aug. 2003)
%

%
%   getting started
      clear all,   close all,   nfig = 0;
%
%   problem data
      alf = 1.0;              % cm^2/s thermal diffusivity
      L = 50;                 % cm     length of bar
      Nx = 101;               % number of x values
      x = linspace(0,L,Nx)';  % vector of points to evaluate function
      nmax = 50;              % max number of nonzero terms
      tol = 0.001;            % tolerance to stop series evaluation
%   
%   calc ln = n*pi/L  and   bn = (4*L/n*pi)*sin(n*pi/2)  for n = 1,3,5,...
      ln = zeros(1,nmax);   bn =  zeros(1,nmax);
      dd = pi/L;    cc = 4*L/(pi^2);    
      for m = 1:nmax
        n = 2*m-1;     ln(m) = n*dd;    cn(m) = cc*sin(n*pi/2)/(n^2);
      end
%
%   now evaluate series expansion for several different times
      tt = [0 25 100 250 500];   Nt = length(tt);   ut = zeros(Nx,Nt);   
%
      for i = 1:Nt
        t = tt(i);   cc = -alf*t;   mrerr = 1.0;   n = 0;   u = zeros(size(x));
        while mrerr > tol   &   n < nmax
          n = n+1;   un = cn(n)*exp(cc*ln(n)*ln(n))*sin(ln(n)*x);   u = u + un;   
          j = find(u);                    % finds indices of nonzero values of u(x)
          mrerr = max(abs(un(j)./u(j)));  % compute max relative error
        end
        ut(:,i) = u;
        disp([' Needed ',num2str(n),' terms for convergence at t = ',num2str(t),' s'])
      end
%
%   plot curves of u(x,t) for various times
       nfig = nfig+1;  figure(nfig)
       v = [0 L 0 25];
       plot(x,ut,'LineWidth',2),axis(v)
       title('HT1D\_SOV1:  Bar Temperature Profile at Various Times (Example 10.1)')
       grid,xlabel('Distance (cm)'),ylabel('Temperature (C)')
       for i = 1:Nt,  gtext(['t = ',num2str(tt(i)),' s']),   end
%
%   end of example