%Temperature distribution in a rod
%See example 9
%Discretization method: Finite difference
% Solution method: SOR
clear all;
a=[];b=[];c=[];d=[];x=[];T=[];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
for k = 1:n % Set coefficients
   a(k)=2.0;
   b(k)=1.0;
   c(k)=1.0;
   d(k)=Qprim*dx^2/lambda;
   T(k) = 30.0; % start value
   x(k) = (k-1)*dx;
end; %
a(1) = 3;
b(1) = 4;
c(1) = 0;
% Solver SOR if omega = 1 Gauss-Siedel
omega=2/(1+sin(pi/nn));
counter = 0;
maxit = 200;
sumres = 1.;
maxres = 1.0e-5;
while ((sumres > maxres)&(counter < maxit) )
   for k = 2:n-1 % SOR
      T(k)=omega*(b(k)*T(k+1)+c(k)*T(k-1)+d(k))/a(k)+(1-omega)*T(k);
   end; % SOR
   d(1) =-T(3);
   T(1)=(b(1)*T(2)+d(1))/a(1);% Insulated end 
   sumres = 0.0
   for k = 2:n-1 % residual
      res=abs(a(k)*T(k)-(b(k)*T(k+1)+c(k)*T(k-1)+d(k)));
      sumres = sumres+res;
   end % residual
   summa = sumres
   counter = counter +1 % counts number of iterations    
end; %while   
%Analytical solution
for k = 1:n
   analytical(k)=Qprim*(L^2-x(k)^2)/(lambda*2)+T(n);
end;
% Plot
hold on;
plot(x,T,'*');
plot(x,analytical,'o');
hold off;
legend('Numerical','Analytical',0);
title('Temperature distribution');