% This program calculates the unsteady heat conduction in a rod
% 
% the thermal diffusivity for aluminum is
%
alpha = 97.1*10-6;  % sq m/sec
L = 4 ; % meters
N = 20; % no of grid elements
delta_x = L/N
delta_t = delta_x^2/alpha
delta_xs = delta_x*delta_x;
itmax = 100;
T0 = 200; % deg C
TL = 100; % deg C
%
% initial temperature profile
for I = 1:N-1
T(I) = TL;
end
%
% calculate matrix elements
a = -alpha*delta_t/(2*delta_x*delta_x);
b = 1 +( alpha*delta_t/(delta_x*delta_x));
c = a;

r1 = [1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 8 8 8 9 9 9];
c1 = [1 2 1 2 3 2 3 4 3 4 5 4 5 6 5 6 7 6 7 8 7 8 9 8 9 10];
v1 = [b c a b c a b c a b c a b c a b c a b c a b c a b c];
r2 = [10 10 10 11 11 11 12 12 12 13 13 13 14 14 14 15 15 15];
c2 = [9 10 11  10 11 12 11 12 13 12 13 14 13 14 15 14 15 16];
v2 = [a  b  c  a  b  c  a  b   c  a  b  c  a  b  c  a  b  c];
r3 = [16 16 16 17 17 17 18 18 18 19 19 ];
c3 = [15 16 17 16 17 18 17 18 19 18 19 ];
v3 = [a  b  c  a  b  c  a  b  c  a  b  ];
r = [r1 r2 r3];
cc = [c1 c2 c3];
v = [v1 v2 v3];
A = sparse(r,cc,v,N-1,N-1);
dd = (1-(alpha*delta_t/delta_xs));
K = 0;
  for I = 1:N+1
   l(I) = (I-1)*delta_x;
  end
l = l';
T = T(:);
TP = [T0;
      T;
      TL];
 if I >= K
   plot(l,TP,'o');
   hold on
   K = K+1;
 end
for I = 1:itmax

 B(1) = -a*T0 +(1-(alpha*delta_t/delta_xs))*T(1) -a*(T0+T(2));
   for J = 2:N-2
      B(J) = dd*T(J) - a*(T(J-1) +T(J+1));
   end
 B(N-1) = -c*TL +dd*T(N-1) -a*(T(N-2) +TL);
 B = B(:);
 T = A\B;
 TP = [T0;
       T;
       TL];
  if I >= K*10
    plot(l,TP);
    hold on
    K = K+1;
  end
 title ('Dynamic Response for Heat Conduction');
 xlabel ('length, m');
 ylabel ('temperature, deg C');
 text(1.25,180,'plot every ten delta time units = 4.15x10-4 sec');
 text(.1,110,'o initial conditions');
end
hold off
