%Temperature distribution in a rod
%See example 9
%Finite volume method
%TDMA
clear all;
T = [];P=[];Q=[];ae=[];aw=[];ap=[];analytical=[];
nn = input('Number of control volumes = ');
n = nn+2;% number of  grid points
L = 0.6;
dx = L/nn; % size of increment
Qprim = 50000; % heat source
lambda = 20; % thermal conductivity
x(1)=0;x(2)=dx/2;
for k = 3:n-1
    x(k)=x(k-1)+dx;
end
x(n)=L;
T(n) = 30;
for k =1:n
   T(k)=30;
end;   
% Solver TDMA
% P(1)=0 and Q(1)=T(1)
for k = 2:n-1
  d(k)=Qprim*dx;
  ae(k)=lambda/dx;
  aw(k)=lambda/dx;
  if (k==2);aw(k)=0;end;
  if (k==n-1);ae(k)=2*lambda/dx;end;
  ap(k)=ae(k)+aw(k);
end;  
P(1) = 0;
Q(1) = T(1);
sumres = 0;
for k = 2:n-1
    P(k) = ae(k)/(ap(k)-aw(k)*P(k-1));
    Q(k) = (aw(k)*Q(k-1)+d(k))/(ap(k)-aw(k)*P(k-1));
end;
for k = n-1:-1:1
    T(k) = P(k)*T(k+1)+Q(k);
end;
T(1)=T(2);
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');