import numpy as np
from matplotlib import pyplot as plt

def f(u):
    rhs = np.zeros(4)
    rhs[0] = u[1]
    rhs[1] = -u[0] + 2.0*u[1]**2/u[0]
    rhs[2] = u[3]
    rhs[3] = (-1.0-2.0*(u[1]/u[0])**2)*u[2] + 4.0*u[1]*u[3]/u[0]
    return rhs

def initialcondition(s):
    u = np.zeros(4)
    u[0] = 1.0/(np.exp(1)+np.exp(-1))
    u[1] = s
    u[2] = 0.0
    u[3] = 1.0
    return u

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

def solvephi(n,s):
    h = 2.0/n
    u = np.zeros((n+1,4))
    u[0] = initialcondition(s)
    for i in range(1,n+1):
        u1 = u[i-1] + 0.5*h*f(u[i-1])
        u[i] = u[i-1] + h*f(u1)
    phi = u[n][0] - 1.0/(exp(1)+exp(-1))
    dphi= u[n][2]
    return phi,dphi,u

n = 100
s = 0.2
p, dp, u = solvephi(n,s)
x = np.linspace(-1.0,1.0,n+1)
xe = np.linspace(-1.0,1.0,1000)
ye = yexact(xe)
plt.plot(x,u[:,0],xe,ye)
plt.legend(("Numerical","Exact"))

n = 100
s = 0.2
maxiter = 100
TOL = 1.0e-8
it = 0
while it < maxiter:
    p, dp, u = solvephi(n,s)
    if np.abs(p) < TOL:
        break
    s = s - p/dp
    it += 1
print "Root = %e" % s
x  = np.linspace(-1.0,1.0,n+1)
xe = np.linspace(-1.0,1.0,1000)
ye = yexact(xe)
plt.plot(x,u[:,0],xe,ye)
plt.legend(("Numerical","Exact"))