function k=atree_bessel_find_range(sigma,p)

# We define temporary function
str = sprintf('function p=atree_bessel_find_range_temp(x)
p=atree_bessel_in_range(%g, min(-x, x), max(-x,x)) - %g; endfunction',
sigma, p);
 eval(str);

# and then solve it for atree_bessel_find_range_temp(x)=0

k0 = [1,0.1,0.01,0.001,sigma];

for i = 1:size(k0)(2)
  [k, info, msg] = fsolve('atree_bessel_find_range_temp',0.1);
  if (info == 1)
    break;
  endif
endfor

k = abs(k);

if (info == 1)
  solution = atree_bessel_in_range(sigma, -k, k);
  if (!atree_compare_round(p, solution))
    fprintf(stderr, "Pr[%g<X<%g]=%g, but Pr[%g<Bessel(%g)<%g]=%g\n", -k, k, p, -k, sigma, k, solution);
    fprintf(stderr, "atree_bessel_find_range didn't get the correct range. fsolve:%s\n",msg);
    perror("fsolve", info);
  endif
else
  fprintf(stderr, "Pr[%g<X<%g]=%g, but Pr[%g<Bessel(%g)<%g]=%g\n", -k, k, p, -k, sigma, k, 
atree_bessel_in_range(sigma, -k,$  fprintf(stderr, "atree_bessel_find_range didn't find the correct range: fsolve 
%s\n",msg);
  perror("fsolve", info);
endif

endfunction