%{
Function that solves a tridiagonal matrix equation
for the backward euler method
%}
function u = tridag(a,b,c,r)
    M = length(a);
    gam = zeros(1,M);
    u = zeros(1,M);
    temp = b(1);
    u(1) = r(1)/temp;
    for j=2:M
        gam(j)=c(j-1)/temp;
        temp = b(j)-a(j)*gam(j);
        if temp==0
            disp('Error, fail');
            return;
        end
        u(j) = (r(j) - a(j)*u(j-1))/temp;
    end
    for j=M-1:-1:1
        u(j) = u(j) - gam(j+1)*u(j+1);
    end
end