%Temperature distribution in a rod
%See example 9
%Finite difference methof
%TDMA
clear all;
T = [];P=[];Q=[];analytical=[];
nn = input('Number of increment = ')
n = nn+1;% number of  grid points
L = 0.6;
dx = L/nn; % size of increment
Qprim = 50000; % heat source
lambda = 20; % thermal conductivity
x = 0:dx:L;
T(n) = 30;
for k =1:n
   T(k)=30;
end;   
% Solver TDMA
% second order appr. dT/dx = 0 gives  3T(1)=4T(2)-T(3)
% to get a tri-diagonal matrix the eqn. T(1) = T(1) is used
% For the results in P(1)=0 and Q(1)=T(1)
a = 2; b = 1; c = 1;
d = Qprim*dx^2/lambda;
counter = 0;
maxit = 200;
sumres = 1.;
maxres = 1.0e-5;
P(1) = 0;
while ((sumres > maxres)&(counter < maxit) )
   Q(1) = T(1);
   sumres = 0;
   for k = 2:n-1
      P(k) = b/(a-c*P(k-1));
      Q(k) = (c*Q(k-1)+d)/(a-c*P(k-1));
      sumres=sumres+abs(a*T(k)-(b*T(k+1)+c*T(k-1)+d));
   end;
   for k = n-1:-1:1
      T(k) = P(k)*T(k+1)+Q(k);
   end;
   T(1)=(4*T(2)-T(3))/3;
   counter=counter +1
end;
for k = 1:n
   analytical(k)=Qprim*(L^2-x(k)^2)/(lambda*2)+T(n);
end;   
hold on;
plot(x,T,'*');
plot(x,analytical,'o');
hold off;
legend('Numerical','Analytical',0);
title('Temperature distribution');