import scipy as N
N.pkgload('special')

import scikits.bvp1lg as bvp

nu = 3.4123
degrees = [2, 1]

def fsub(x, z):
    """The equations"""
    u, du, v = z     # it's neat to name the variables
    return N.array([-du/x + (nu**2/x**2 - 1)*u, x**(nu+1) * u])

def gsub(z):
    """The boundary conditions"""
    u, du, v = z
    return N.array([u[0] - N.special.jv(nu,   1),
                    v[1] - 5**(nu+1) * N.special.jv(nu+1, 5),
                    u[2] - N.special.jv(nu,   10)])

boundary_points = [1, 5, 10]
tol = [1e-5, 0, 1e-5]

solution = bvp.colnew.solve(
    boundary_points, degrees, fsub, gsub,
    is_linear=True, tolerances=tol,
    vectorized=True, maximum_mesh_size=300)

# Analytical solution is of course known here, so check it

x = solution.mesh
assert N.allclose(solution(x)[:,0], N.special.jv(nu, x),
                  rtol=2e-3, atol=1e-8)
assert N.allclose(solution(x)[:,2], x**(nu+1)*N.special.jv(nu+1, x),
                  rtol=2e-3, atol=1e-8)

# Plot it

import pylab

pylab.plot(x, N.special.jv(nu, x), 'g.',
           x, solution(x)[:,0], 'b-')
pylab.title('$u$')
pylab.savefig('bessel-solution.png')
pylab.show()

pylab.plot(x, N.special.jv(nu+1, x)*x**(nu+1), 'g.',
           x, solution(x)[:,2], 'b-',)
pylab.title('$v$')
pylab.savefig('bessel-solution-v.png')
pylab.show()