import scipy as N
N.pkgload('special')
import scipy as N
import scikits.bvp1lg as bvp

C = None

# The exact solution is known

def exact_solution(x):
    return N.arccosh(N.sqrt((1.+C)/C)*N.cosh(N.sqrt(C)*(x - .5)))

def f(x, z):
    """Right-hand side"""
    u, du = z
    return N.array([N.cosh(u)/N.sinh(u)**3])

def df(x, z):
    """Jacobian of f"""
    u, du = z
    return N.array([[-(1 + 2*N.cosh(u)**2)/N.sinh(u)**4, 0*x]])

def g(z):
    """Boundary conditions"""
    u, du = z
    v0 = exact_solution(0)
    v1 = exact_solution(1)
    return N.array([u[0] - v0,
                    u[1] - v1])

def dg(z):
    """Jacobian of g"""
    return N.array([[1, 0],
                    [1, 0]])

def guess(x):
    z = N.zeros((2,) + N.asarray(x).shape)
    dm = N.zeros((1,) + N.asarray(x).shape)
    z[0] = 1
    return z, dm

# Solve the problem for multiple values of C using simple continuation.

Cvals = N.r_[N.linspace(1, 1.69, 5),
             N.linspace(1.69, 1.72, 10), # branch point here?
             N.linspace(1.8, 500, 40)]

x = N.linspace(0, 1, 501)

solution = guess
for Cval in Cvals:
    C = Cval
    print C
    solution = bvp.colnew.solve([0, 1], [2],
                                f, g,
                                dfsub=df, dgsub=dg,
                                initial_guess=solution,
                                tolerances=[1e-5, 1e-5],
                                maximum_mesh_size=1500,
                                vectorized=True,
                                coarsen_initial_guess_mesh=True)
    assert N.allclose(exact_solution(x), solution(x)[:,0], rtol=1e-5)

# Plot it

import pylab

x = solution.mesh

pylab.plot(x, exact_solution(x), 'g.',
           x, solution(x)[:,0], 'b-')
pylab.title('$u$')
pylab.savefig('nonlin-solution.png')
pylab.show()