I have tried boost::math::cyl_bessel_k(x,y) * exp(y). In most cases, this is equal to Matlab's scaled besselk(x,y,1). But in some cases (e.g., x=1, y=2000) when both besselk(x,y)=0 and boost::math::cyl_bessel_k(x,y)=0, Matlab's scaled version besselk(x,y,1) gives me different values varies around 10^-3. But boost::math::cyl_bessel_k(x,y) * exp(y) returns -nan. I'd like to find an equivalent statement for Matlab's besselk(x,y,1). How can I handle this?
I'm not seeing anything in Boost that does what you need (though you might be able to implement it yourself by using lower-level functions). As you've found out, the scaled Bessel functions are not computed simply by multiplying exp(z). The GSL appears to have incorporated this functionality, e.g., gsl_sf_bessel_Knu_scaled. For an "exact equivalent," you might look at the paper and code by Amos, e.g., CBESK. Both Matlab and Octave appear to use this implementation. However, the code is written in Fortran, so you'd need to translate it or put a wrapper around it (this project appears to have done that so it might be useful – there are others out there too).
You may also be able to use Matlab's Coder and codegen to output something as well.