import numpy

def backsubstitute(A, b):
    n = len(b)
    i,j = A.shape
    if i != n or j != n:
        raise ValueError("A should be shape (%d, %d), but is %s" % (n,n,A.shape))

    x = b.copy()

    for i in range(n-1,-1,-1):
        x[i] = b[i]
        for j in range(n-1, i, -1):
            x[i] -= A[i,j] * x[j]
        x[i] /= A[i,i]

    return x

def gauss_eliminate(A_in, b_in):
    # Kopien der Matrix und des Vektors anlegen
    A = A_in.copy()
    b = b_in.copy()

    # Gausselimination durchfuehren
    # ....

    return A, b

def solve(A, b):
    Ares, bres = gauss_eliminate(A, b)
    x = backsubstitute(Ares, bres)
    return x

# Sample solutions
A1 = array([[ 1.        ,  0.2       ,  0.04      ,  0.008     ],
       [ 1.        ,  0.53333333,  0.28444444,  0.1517037 ],
       [ 1.        ,  0.86666667,  0.75111111,  0.65096296],
       [ 1.        ,  1.2       ,  1.44      ,  1.728     ]])
x1 = array([ 1.27785932, -0.17219435, -8.75422788,  7.22564732])
b1 = array([ 0.95105652, -0.20791169, -0.74314483,  0.95105652])

A2 = array([[ 0.        ,  0.2       ,  0.04      ,  0.008     ],
       [ 1.        ,  0.53333333,  0.28444444,  0.1517037 ],
       [ 1.        ,  0.86666667,  0.75111111,  0.65096296],
       [ 1.        ,  1.2       ,  1.44      ,  1.728     ]])
x2 = array([ -0.85418404,   8.06213979, -18.74818115,  11.0694755 ])
b2 = array([ 0.95105652, -0.20791169, -0.74314483,  0.95105652])