def finite_dif(p, q, r, a, b, alpha, beta, N):
'''
Approximate the solution to the boundary value problem
y'' = p(x)y' + q(x)y + r(x), r<=a<=b, y(a)=alpha, y(b)=beta
INPUT
Functions p(x), q(x), r(x); limit conditions a, b; boundary conditions alpha, beta; integer N.
OUTPUT
Approximations w[i] of y[xi] for each i = 0, 1, ..., N+1.
'''
an = []
bn = []
cn = []
dn = []
h = (b-a)/(N+1)
x = a+h
an.append(2 + (h**2)*q(x))
bn.append(-1+(h/2.0)*p(x))
cn.append(0)
dn.append(-(h**2)*r(x)+(1+(h/2.0)*p(x))*alpha)
for i in range(2, N):
x = a+i*h
an.append(2+(h**2)*q(x))
bn.append(-1+(h/2)*p(x))
cn.append(-1-(h/2)*p(x))
dn.append(-(h**2)*r(x))
x = b-h
an.append(2+(h**2)*q(x))
bn.append(0)
cn.append(-1-(h/2)*p(x))
dn.append(-(h**2)*r(x)+(1-(h/2.0)*p(x))*beta)
A = numpy.zeros([N,N])
A[0][0] = an[0]
A[0][1] = bn[0]
for j in range(1,N-1):
A[j][j-1] = cn[j]
A[j][j] = an[j]
A[j][j+1] = bn[j]
A[N-1][N-2] = cn[N-1]
A[N-1][N-1] = an[N-1]
B = matrix(A)
d = vector(dn)
wn = B.solve_right(d)
xn = []
for k in range(1,N+1):
xn.append(a+k*h)
print '%6s %6s'%("X", "Wi")
print '%6s %6s'%(a, alpha)
for l in range(len(xn)):
print '%6s %6s'%(xn[l],wn[l])
print '%6s %6s'%(b, beta)
return