Include a function evaluation too:
def f(x): return x**2 + 1 # Scalar version for i in range(1, n-1): u2[i] = u[i-1] - 2*u[i] + u[i+1] + f(x[i]) # Vectorized version u2[1:-1] = u[:-2] - 2*u[1:-1] + u[2:] + f(x[1:-1])