[jprimemk,JofJ]=BessDerivZerosBisect(mMAX,kMAX)

function [jprimemk,JofJ]=BessDerivZerosBisect(mMAX,kMAX)
%This routine calculates a table of the (k) zeros of the Derivative of
%the m-th order Bessel function J'm.
%   m , k  =  1 , 2 , 3 , ...
%Use the bisection algorithm, because it gives the desired roots.
%Estimates for bracketing intervals are taken from the
%asymptotic forms for the roots, in Abramowitz and Stegun.
%
%jprimemk  are the roots, and  JofJ  is Bessel function evaluated there.
%
%Routine written by Larry Forbes, University of Tasmania.
%
jprimemk=zeros(mMAX,kMAX);
JofJ=zeros(mMAX,kMAX);
options=[];
for m=1:mMAX
    for k=1:kMAX
        %give the asymptotic guess at the kth zero of J'm .
        if k == 1
            oneterm=m+0.8086165*m^(1/3)+0.072490*m^(-1/3)-0.05097*m^(-1);
            oneterm=oneterm+0.0094*m^(-5/3);
            asymptroot=oneterm;
        else
            betapr=(m/2+k-3/4)*pi;
            mu=4*m^2;
            asymptroot=betapr-(mu+3)/(8*betapr);
            asymptroot=asymptroot-4*(7*mu^2+82*mu-9)/(3*(8*betapr)^3);
        end
        %bracket the root, and check that the root is bracketed.
        errmess='original interval does not bracket root';
        if k == 2
            Aleft=asymptroot-pi/1.2;
        elseif k == 3
            Aleft=asymptroot-pi/2;
        else
            Aleft=asymptroot-pi/(3+0.02*k);
        end
        Aright=asymptroot+pi/3;
        JprimeL=besselderiv(Aleft,m);
        JprimeR=besselderiv(Aright,m);
        if JprimeL*JprimeR > 0
            disp(errmess)
            m
            k
        end
        Amiddle=(Aleft+Aright)/2;
        JprimeM=besselderiv(Amiddle,m);
        while abs(JprimeM) > 1.e-6
            if JprimeL*JprimeM < 0
                Aright=Amiddle;
                JprimeR=besselderiv(Aright,m);
                Amiddle=(Aleft+Aright)/2;
                JprimeM=besselderiv(Amiddle,m);
            else
                Aleft=Amiddle;
                JprimeL=besselderiv(Aleft,m);
                Amiddle=(Aleft+Aright)/2;
                JprimeM=besselderiv(Amiddle,m);
            end
        end
        jprimemk(m,k)=Amiddle;
        JofJ(m,k)=JprimeM;
    end
end
function derivvalue=besselderiv(x,m)
thederiv=(BesselJ(m-1,x)-BesselJ(m+1,x))/2;
derivvalue=thederiv;