# 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:
# [a0 d0  0  0  0
#  c1 a1 d1  0  0
#  0  c2 a2 d2  0
#  0  0  c3 a3 d3
#  0  0  0  c4 a4]

import numpy as np

def thomas(a,d,c,b) :
    n = len(a)
    x = np.zeros(n)
    for i in range(1,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]

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

    return x