import numpy as np
from matplotlib import pyplot as plt

def yexact(x):
    return 1.0/(np.exp(x)+np.exp(-x))

def f(y,dy):
    return -y + 2.0*dy**2/y

def phi(y):
    n   = len(y) - 1
    h   = 2.0/n
    res = np.zeros(n-1)
    for i in range(1,n):
        dy = (y[i+1] - y[i-1])/(2.0*h)
        res[i-1] = (y[i-1]-2.0*y[i]+y[i+1])/h**2 - f(y[i],dy)
    return res
        
def dphi(y):
    n   = len(y) - 1
    h   = 2.0/n
    res = np.zeros((n-1,n-1))
    for i in range(1,n):
        dy = (y[i+1] - y[i-1])/(2.0*h)
        if i > 1:
            res[i-1,i-2] = 1.0/h**2 + 4.0*dy/y[i]
        res[i-1,i-1] =-2.0/h**2 + 1.0 - 2.0*(dy/y[i])**2
        if i < n-1:
            res[i-1,i  ] = 1.0/h**2 - 4.0*dy/y[i]
    return res

n = 100

# Initial guess for y
y = np.zeros(n+1)
y[:] = 1.0/(exp(1) + exp(-1))

maxiter = 100
TOL = 1.0e-6
it = 0
while it < maxiter:
    b = phi(y)
    if np.linalg.norm(b) < TOL:
        break
    A = dphi(y)
    v = np.linalg.solve(A,b)
    y[1:n] = y[1:n] - v
    it += 1
print "Number of iterations = %d" % it
x = np.linspace(-1.0,1.0,n+1)
plt.plot(x,y,'o',x,yexact(x))
plt.legend(("Numerical","Exact"));