Calculate second order bessel functions in java with non-integer order

Question

I am trying to calculate second order bessel functions (Y_v(x)) using Java. Where v, the order of the bessel function, is non integer.

Apache, I know, has a function to calculated first order bessel functions (J_v(x))

import org.apache.commons.math3.special.BesselJ;

BesselJ bj = new BesselJ(10.5);
double bjr = bj.value(r);

Which can calculate non-integer first order bessel functions (but not negative order functions). I do not believe it has an equivalent for second order.

I have found a formula for calculating second order bessel functions (equation 6.4.2 of Numerical Recipies among other places) which calculates second order bessel functions from first order functions.

However this equation requires the first order bessel functions to accept negative v values which the apache BesselJ function does not.

I have written a program to calculate BesselJ functions whch accept negative values using the formula (Numerical Recipes 6.4.1)

import org.apache.commons.math3.special.Gamma;
import org.apache.commons.math3.util.ArithmeticUtils;

private static int lmax;
private static double[][] Gammas;
private static int besselcount;

public static Double Bessel_J(double m, double x){

    double BJ = 0.0;
    for (int k = 0; k < besselcount;k++){
        double kk = (double) k;     
        double llmm = (double) lmax;
        double A = Math.pow(x/2.,(2.*kk)+m);
        double B = Gammas[(int) (m+llmm+0.5)][k];
        if (B==0.0||B==-0.0){break;}
            BJ = BJ + A*B;
        }

        return BJ;
}

public static void main(String[] args) throws Exception {

    besselcount = 80;
    lmax = 20;

    Gammas = new double[2*lmax+2][besselcount+1];
    for (int n = 0;n <= ((2*lmax)+1);n++){
        for (int m = 0;m <=besselcount;m++) {
            double nn = (double) n;
            double mm = (double) m;
            double llmm = (double) lmax;
            Gammas[n][m] =   Math.pow(-1,m)/(ArithmeticUtils.factorialDouble(m)*Gamma.gamma(1.+nn-llmm-0.5+mm));
        }
    }

    for (double x = 0.02; x < 50.0 ; x = x + 0.02){
        double bj = Bessel_J(-10.5, x);
        System.out.println(x+" "+bj);
    }

}

(Having the gammas precalculated and stored in an array is just because of how the bessel function is used in the rest of the code, to avoid recalculating this unnecessarily)

The function I have written works for low x values but then looses stability. For the v = -10.5 value I have used after around x = 30. As shown in the figure.

On the wolfram alpha website is says mathematica uses the same equations to calculate bessel functions and mathematica's BesselJ function can calculate

Show[Plot[BesselJ[-10.5, x], {x, 0, 50}], ImageSize -> Full]

Without loosing stability.

Any advice on either how to correct my functions or alternative ways to calculate second order bessel functions using Java would be very much appreciated.

Please say if any clarification would be useful.

Thank you