from scipy import optimize, special
from numpy import *
from pylab import *
 
x = arange(0, 10, 0.01)
 
for k in arange(0.5, 5.5):
    y = special.jv(k, x)
    plot(x, y)
    f = lambda x:-special.jv(k, x)
    x_max = optimize.fminbound(f, 0, 6)
    plot([x_max], [special.jv(k, x_max)], 'ro')

title('Different Bessel functions and their local maxima')
show()
