# function thomas
# Inputs are vectors a,d,c,b where a,d, and c are the diagonal, upper diagonal, and
#    lower diagonals of A, and b is the rhs vector:
# [a1 d1  0  0  0
#  c2 a2 d1  0  0
#  0  c3 a3 d3  0
#  0  0  c4 a4 d4
#  0  0  0  c5 a5]

function thomas(a,d,c,b)
    n = length(a)
    x = zeros(n)
    for i in 2:n
        a[i] = a[i] - c[i]/a[i-1]*d[i-1]
        b[i] = b[i] - c[i]/a[i-1]*b[i-1]
    end

    x[end] = b[end]/a[end]
    for i in n-1:-1:1
        x[i] = (b[i]-d[i]*x[i+1])/a[i]
    end

    return x
end