#%%
from numpy import array, empty, linalg, zeros
from matplotlib import pyplot
from scipy import interpolate

def Miguel_S5_Aufg2(x, y, xx):
    n = len(x)
    a = empty((n - 1, 1))
    b = empty((n - 1, 1))
    c = empty((n + 1, 1))
    d = empty((n - 1, 1))
    h = empty((n - 1, 1))

    for i in range(0, n - 1):
        a[i] = y[i]
        h[i] = x[i + 1] - x[i]

    c[0] = 0
    c[n] = 0

    A = zeros((n - 2, n - 2))
    z = empty((n - 2, 1))

    for i in range(0, n - 2):
        if i == 0:
            A[i][0] = 2 * (h[0] + h[1])
            A[i][1] = h[1]
            z[i] = 3 * (((y[2] - y[1]) / h[1]) - ((y[1] - y[0]) / h[0]))
        elif i == n - 3:
            A[i][i - 1] = h[i - 1]
            A[i][i] = h[i - 1] + h[i]
            z[i] = 3 * (((y[i + 1] - y[i]) / h[i]) - ((y[i] - y[i - 1]) / h[i - 1]))
        else:
            A[i][i - 1] = h[i - 1]
            A[i][i] = 2 * (h[i - 1] + h[i])
            A[i][i + 1] = h[i]
            z[i] = 3 * (((y[i + 1] - y[i]) / h[i]) - ((y[i] - y[i - 1]) / h[i - 1]))

    c[1:n - 1,:] = linalg.solve(A, z)	

    for i in range(0, n - 1):
        b[i] = ((y[i + 1] - y[i]) / h[i]) - ((h[i] / 3) * (c[i + 1] + 2 * c[i]))
        d[i] = (1 / (3 * h[i])) * (c[i + 1] - c[i])

    yy = empty(len(xx))

    for i in range(0, len(yy)):
        j = 0
        value = xx[i]

        while x[j] < value:
            j += 1

        yy[i] = y[j - 1] + b[j - 1] * (value - x[j - 1]) + c[j - 1] * (value - x[j - 1]) ** 2 + d[j - 1] * (value - x[j - 1]) ** 3

    pyplot.scatter(xx, yy)
    return yy

# %%

x = array([4, 6, 8, 10])
y = array([6, 3, 9, 0])
xx = array([4.5, 6.5, 9.0])

Miguel_S5_Aufg2(x, y, xx)
pyplot.scatter(xx, interpolate.CubicSpline(x, y)(xx))

# %%