Algebra_solver with integrate_1d containing modified_bessel_first_kind

General
,
1

I created a function for the PDF of an uncorrelated bivariate normal distribution as a function of radius (equation 17 from A Toolbox for the Radial and Angular Marginalization of Bivariate Normal Distributions (arxiv.org)). My goal is to integrate that, and solve to return the quantile associated with a given cumulative probability. For example, in the code below I want to find the radius r that contains p = 50\% of the data given an uncorrelated bivariate normal distribution centered at the origin with standard deviations \sigma_x = 1, and \sigma_y = 1.

I have tested this and can integrate with no problem, but trying to solve for p is a problem that seems to be due to modified_bessel_first_kind. If I comment out the line with modified_bessel_first_kind it works, but otherwise I get the error: Exception: gradient_of_f: The gradient of f is nan for parameter 0. Is there any way that I might be able to work around this?

functions{
        real r_marginal_polar_bivariate_normal_pdf(
          real x, real xc, real[] theta, real[] x_r, int[] x_i) {
                    
          // define my parameters
          real sx = theta[1];
          real sy = theta[2];
          real prob;

          prob = (x/(sx*sy)) *
                 exp(-(sx^2 + sy^2)/((2*sx*sy)^2) * x^2) *
                 modified_bessel_first_kind(0, -(sx^2 - sy^2)/((2*sx*sy)^2) * x^2 );
          return prob;
        }

        vector system(vector y, vector theta, real[] x_r, int[] x_i){
          vector[1] root;
          real sx;
          real sy;
          real p;
          sx = theta[1];
          sy = theta[2];
          p = theta[3];
          root[1] = integrate_1d(r_marginal_polar_bivariate_normal_pdf, 0, y[1], {sx, sy},
            rep_array(0.0, 0), rep_array(0, 0), 1e-8)-p;
          return root;
        }       
}
generated quantities {
  real x_r[0];
  int x_i[0];

  vector[1] y;
  vector[1] y_guess;
  vector[3] theta;
  theta[1] = 1;       // sx
  theta[2] = 1;       // sy
  theta[3] = .5;      // p
  y_guess[1] = 1.17;
  
  // find the quantile of r corresponding with p
  y = algebra_solver(system, y_guess, theta, x_r, x_i, 1e-8, 1e-5, 1e4);
}