Overflow of Bessel functions in numerical solution of characteristic equation with double-precision for large-core and high numerical aperture optical fibers is circumvented by evaluating the ratio of the nth order to the (n-1)th order modified Bessel functions of the second kind using a small argument approximation and by evaluating the ratio of nth order to the (n-1)th order Bessel function of the first kind using Lentz-Thompson method. The evaluation error of the two methods is not larger than the minimum value of the floating-point system used. Characteristic equations of multimode silica optical fiber and large core polymer optical fiber are solved using this algorithm and compared with commercial optical software (IFO-Grating) from Optiwave Corp. The results for the silica optical fiber are good agreement with the commercial software. Authors' algorithm has solved all guided modes for the large core optical fiber and results for the guided modes with angular mode number smaller than 70 are same to that of IFO-Grating, but IFO-Grating failed to solve the eigenvalue for guided modes with angular mode number larger than 70.