%{
  Solution of tridiagonal linear systems with Thomas' algorithm

 clear all;
 l = [0;-1;-1]; d = [2;2;2]; u = [-1;-1;0]; b = [1;2;3];
 A = spdiags([[l(2:end);nan],d,[nan;u(1:end-1)]],[-1,0,1],3,3);
 xm = A\b;
 xt = thomas(l,d,u,b);
%}
function x = thomas(l,d,u,b)

N = numel(d);

% Gaussian elimination (foreward)
for j=1:N-1
    m = l(j+1)/d(j);
    d(j+1) = d(j+1)-m*u(j);
    b(j+1) = b(j+1)-m*b(j);
end


% Backward substitution
x(N,1) = b(N)/d(N);
for j=N-1:-1:1
    x(j) = (b(j) - u(j)*x(j+1))/d(j);
end

end