import numpy as np
import matplotlib.pyplot as plt

n = 40 #Number of interior points
h = 1./float(n+1) #size of partition
x = np.linspace(h,1,n+1) #x[i] = (i+1)h

#Linear system to be solved: Av = b, A is (n+1 x n+1). Index runs from 0,...,n
#Define A_1
A_1 = np.zeros((n+1,n+1))
for i in range(0,n):
    A_1[i][i] = 2
A_1[n][n] = 1
for i in range(0,n):
    A_1[i+1][i] = -1
    A_1[i][i+1] = -1         

#Define A_2
A_2 = np.zeros((n+1,n+1))
for i in range(0,n+1):
    A_2[i][i] = 2
for i in range(0,n):
    A_2[i][i+1] = -1
    A_2[i+1][i] = -1
A_2[n][n-1] = -2

#Define b_1
b_1 = -(h**2)*np.exp(x-1)
b_1[n] = h

#Define b_2
b_2 = -(h**2)*np.exp(x-1)
b_2[n] = b_2[n] + 2*h

#Solving system
v_1 = np.linalg.solve(A_1,b_1)
v_2 = np.linalg.solve(A_2,b_2)

#Computing error
u = np.exp(-1)*(np.exp(x)-1)
E_1 = np.amax(np.fabs(v_1-u))
E_2 = np.amax(np.fabs(v_2-u))
print 'E_1: ', E_1/h
print 'E_2: ', E_2/h**2

#Plot numerical solution
plt.plot(x,v_1)
plt.plot(x,v_2)

#Exact solution along a finer grid
x = np.linspace(0,1,100)
u = np.exp(-1)*(np.exp(x)-1)
plt.plot(x,u, '--')

plt.show()