import numpy as np
import numdifftools as nd
import math
import matplotlib.pyplot as plt

#enter h*i where there is a t
def f1(h,i):
    return -1
def f2(h,i):
    return ((h**2)*(math.pi**2)/9)-2
def f3(h,i):
    return -1

#if there is a number outside f then put in z[i]
#interval [a,b], boundary values[c,d], step size h
def g(a, b, c, d, n):
    #The Matrix
    h = float(abs(b - a)/(n + 1))
    w = np.zeros(shape = (n, n))
    i = 0
    w[i, i] = f2(h,i)
    w[i,i+1] = f3(h,i)
    i = 1
    while i < n-1:
        w[i, i - 1]= f1(h,i)
        w[i, i] = f2(h,i)
        w[i, i + 1] = f3(h,i)
        i = i + 1
        pass
    w[i, i - 1]= f1(h,i)
    w[i, i] = f2(h,i)

    #solution vector
    z = np.zeros(shape = (n, 1))
    z[0] = c
    i = 1
    while i<n-1:
        z[i]=0
        i = i+1
        pass
    z[n-1] = d
    p=np.linalg.solve(w,z)

    #finished vector
    y = np.zeros(shape = (n+2, 1))
    y[0] = c
    y[n+1] = d
    i=0
    while i<n:
        y[i+1]=p[i]
        i = i+1
        pass
    return y

#Exact Solution
def exact(t):
    return 3*math.sin(t*math.pi/3)-math.cos(t*math.pi)/3
def actual(a, b, c, d, n):
    h = float(abs(b - a)/(n + 1))
    w = np.zeros(shape = (n+2, 1))
    w[0] = c
    w[n+1] = d
    i = 1
    while i < n+2:
        w[i] = exact(h*i)
        i = i + 1
        pass
    return w

#Answers
e39 = g(0.0,1.5,-1.0,3.0,39)
a39 = actual(0.0,1.5,-1.0,3.0,39)
print(e39)
print(a39)
plt.plot(e39)
plt.plot(a39)
print(a39-e39)

e19 = g(0.0,1.5,-1.0,3.0,19)
a19 = actual(0.0,1.5,-1.0,3.0,19)
print(e19)
print(a19)
plt.plot(e19)
plt.plot(a19)
print(a19-e19)

e9 = g(0.0,1.5,-1.0,3.0,9)
a9 = actual(0.0,1.5,-1.0,3.0,9)
print(e9)
print(a9)
plt.plot(e9)
plt.plot(a9)
print(a9-e9)

plt.show()