Bessel functions of the first kind, \(J_l\) where \(l\) is an integer,
are found in `scipy.special.jn`. This is a function that takes two arguments:
the first one is the value of \(l\) and the second one is the value of the
point at which to evaluate the functions. If one is interested in \(J_0\)
or \(J_1\) only, then there are faster implementations of these functions
as `j0` and `j1`.

If the goal is to find the zeros of the Bessel function, then rather than using
a root finding routine (such as `scipy.optimize.newton`) it is more efficient
to call `scipy.special.jn_zeros`. This function also takes two arguments, the
first one is the value of \(l\) and the second one is the number of zeros
to display (it starts at the first positive zero).

For Bessel functions of half-integer order, so called *spherical Bessel
function*, the situation is a little bit more complicated. The appropriate
function is `scipy.special.sph_jn`. To compute \(J_{l + 1/2}\), the first
argument to `sph_jn` is the value of \(l\). The second argument is the
value of the point at which to evaluate the function. Where it gets complicated
is in the return value of the function. `sph_jn` returns a pair of arrays.
The first array contains the values of the function, and the second array
contains the values of the derivative. Each array has length \(l + 1\) and
contains a value for all spherical Bessel functions with \(l\) smaller than
or equal to the supplied value.

It is useful to wrap this function into something friendlier:

def j3o2(x): """Compute the spherical Bessel function J_{3/2}.""" return sph_jn(1, x)[0][1] def j5o2(x): """Compute the spherical Bessel function J_{5/2}.""" return sph_jn(2, x)[0][2]

To obtain the zeros of the spherical Bessel function, one has to call a root
finding procedure. For instance, `scipy.optimize.newton` does the job quite
well.

John D. Cook has a great page about Bessel functions and another one about how to work with the implementation in Scipy.