title: Discrepancy between analytic solution and solution by relaxation method
body: <p>So I am trying to solve the differential equation $\frac{d^2y}{dx^2} = -y(x)$ subject to boundary conditions y(0) = 0 and y(1) = 1 ,the analytic solution is y(x) = sin(x)/sin(1).</p>
  <p>I am using three point stencil to approximate the double derivative.
  The curves obtained through these ways should match at least at the boundaries ,but my solutions have small differences even at the boundaries.
  I am attaching the code, Please tell me what is wrong.</p>
  <pre><code>import numpy as np
  import scipy.linalg as lg
  from scipy.sparse.linalg import eigs
  from scipy.sparse.linalg import inv
  from scipy import sparse
  import matplotlib.pyplot as plt
  
  a = 0
  b = 1
  N = 1000
  h = (b-a)/N
  r = np.arange(a,b+h,h)
  
  y_a = 0
  y_b = 1
  
  def lap_three(r):
      h = r[1]-r[0]
      n = len(r)
     
      
      M_d = -2*np.ones(n)
     #M_d = M_d + B_d
      O_d = np.ones(n-1)
      mat = sparse.diags([M_d,O_d,O_d],offsets=(0,+1,-1))
     #print(mat)
      return mat
  
  
  def f(r):
      h = r[1]-r[0]
      n = len(r)
      return -1*np.ones(len(r))*(h**2)
  
  def R_mat(f,r):
      r_d = f(r)
      R_mat = sparse.diags([r_d],offsets=[0])
     #print(R_mat)
      return R_mat
  
  #def R_mat(r):
   #   M_d = -1*np.ones(len(r))
      
      
  def make_mat(r):
      main = lap_three(r) - R_mat(f,r)
      return main
  
  main = make_mat(r)
  main_mat = main.toarray()
  print(main_mat)
  '''
  eig_val , eig_vec = eigs(main, k = 20,which = 'SM')
  
  #print(eig_val)
  Val = eig_vec.T
  plt.plot(r,Val[0])
  '''
  main_inv = inv(main)
  inv_mat = main_inv.toarray()
  #print(inv_mat)
  #print(np.dot(main_mat,inv_mat))
  n = len(r)
  B_d = np.zeros(n)
  B_d[0] = 0
  B_d[-1] = 1
  #print(B_d)
  #from scipy.sparse.linalg import spsolve
  A = np.abs(np.dot(inv_mat,B_d))
  plt.plot(r[0:10],A[0:10],label='calculated solution')
  
  real = np.sin(r)/np.sin(1)
  
  plt.plot(r[0:10],real[0:10],label='analytic solution')
  plt.legend()
  #plt.plot(r,real)
  #plt.plot(r,A)
  
  '''diff = A-real
  plt.plot(r,diff)'''
  </code></pre>