I am attempting to solve the following ODE problem:

$$-u''+ u = x$$ $$u(0) = 0$$ $$u'(1) = -u(1)$$

The exact solution is: $u(x) = e^{-x-1} - e^{x-1} + x$

I have a Dirichlet at $x = 0$ and a Robin condition at $x = 1$. I refer to the exact solution with "little" $u$, and the approximate solution with "big" $U_i$, where the subscript $i$ refers to "point number": if I have $m$ sub-intervals, then I will have $m+1$ points indexed from $0$ to $m$.

In order to implement the Dirichlet condition, I eliminated $u_0$ from the finite difference matrix, and modified the right hand side, but in this case I didn't really modify the RHS since $u(0) = 0$, thus $U_0 = 0$.

$$u''(0+h) \approx \frac{U_0 - 2U_1 + U_2}{h^2} = \frac{- 2U_1 + U_2}{h^2}$$

In order to implement the Robin condition, I decided to use the "ghost point" method along with a centred difference approximation of the first derivative ($U_m = u(1)$, if there are $m$ points).

$$u'(1) = -u(1),\;\therefore u'(1) \approx \frac{U_{m+1} - U_{m-1}}{2h} = -U_m$$ $$U_{m+1} = U_{m-1} - 2hU_{m}$$ $$\therefore u''(1) \approx \frac{U_{m-1} - 2U_m + U_{m+1}}{h^2} = \frac{2U_{m-1} + (-2 - 2h)U_m}{h^2}$$

I then implemented the above approximation using some simple Python code (you should be able to run it too), and I have included comments liberally to aid in readability:

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt

#Want to solve ODE: -u'' + u = x, u(0) = 0, u'(1) = -u(1)
#Actual solution: exp(-X-1)
# a and b mark the interval over which we want to solve our ODE
a = 0
b = 1
num_subintervals = 4
num_points = num_subintervals+1

#size of each subinterval
h = (b - a)/num_subintervals

X = np.linspace(a, b, num_points)

def f(x):
    return x
vectorized_f = np.vectorize(f)

# actual solution: -u'' + u = x, u(0) = 0, u'(1) = -u(1)
u = np.exp(-X-1) - np.exp(X-1) + X

#Implementing Dirichlet condition at u(0)
F = vectorized_f(X)    
u_0 = 0
F[1] = F[1] - (u_0/(h**2))
#We have eliminated U_0 from our equations due to implementation of BCs
F = F[1:]

#eliminated U_0, but U_n is still in the equation, thus num_points-1 to solve for
N = num_points-1
# -2 along main diagonal, 1 along +1 and -1 diagonal
D2 = np.diag((-2/(h**2))*np.ones(N), 0) + np.diag((1/(h**2))*np.ones(N-1), 1) + np.diag((1/(h**2))*np.ones(N-1), -1)

#D1(U_n) = (U_n+1 - U_n-1)/2h = -U_1
# -2h*U_n + U_n-1 = U_n+1
# D2(U_n) = (U_n-1 - 2*U_n + U_n+1)/h**2 = (U_n-1 - 2*U_n - 2*h*U_n + U_n-1)/h**2
# D2(U_n) = (-2*U_n - 2*h*U_n + 2*U_n-1)/h**2
D2[N-1, N-2] = 2/(h**2) #Modifying U_n-1 coefficient (index N-2) in equation for U''_n
D2[N-1, N-1] = (-2 - 2*h)/(h**2) #Modifying U_n (index N-1) in equation for U''_n

I = np.ones((num_points-1, num_points-1))

#LHS
A = -1*D2 + I

U = np.concatenate(([u_0], np.linalg.solve(A, F)))

error = np.linalg.norm(U - u, ord=np.inf)

print "error: ", error
plt.plot(X, U, label="approximate")
plt.plot(X, u, label="exact")
plt.legend(loc="best")
plt.show()

I get the following image:

And the error does not decrease as I increase number of sub-intervals. So, clearly, I have an error, but I am not sure where. Could you help?