% Heated Rod Example
% Coded by Nigel F. Reuel on 10.17.17
% This code solves the heated rod question with '\' and tridiag algorithm
% See visual slides for derivation of this coefficient matrix
A = [-2.04 1 0 0; 1 -2.04 1 0; 0 1 -2.04 1; 0 0 1 -2.04];
b = [-40.8; -.8; -.8; -200.8];
tic 
Tvec = A\b;
time1 = toc
% Now solve this problem with Tridiag
e = [0; 1; 1; 1];
f = ones(1,4)*-2.04;
g = [1; 1; 1; 0];
r = b;
tic
Tvec2 = Tridiag(e,f,g,r);
time2 = toc
% Plot temp profile
T = [40 Tvec2 200];
X = 0:2:10;
plot(X,T)
% NOTE tridiag is not much faster for this TINY matrix, but when you have a
% real finite element system with MANY interior nodes this simplified gauss
% elimination approach is VERY powerful.

function x = Tridiag(e,f,g,r)
% Tridiag: Tridiagonal equation solver banded system
% x = Tridiag(e,f,g,r): Tridiagonal system solver.
% input:
% e = subdiagonal vector
% f = diagonal vector
% g = superdiagonal vector
% r = right hand side vector
% output:
% x = solution vector
n=length(f);
% forward elimination
for k = 2:n
factor = e(k)/f(k-1);
f(k) = f(k) - factor*g(k-1);
r(k) = r(k) - factor*r(k-1);
end
% back substitution
x(n) = r(n)/f(n);
for k = n-1:-1:1
x(k) = (r(k)-g(k)*x(k+1))/f(k);
end
end