We present an algorithm to evaluate the modified Bessel functions I v and K v of integral and half-integral order based onthe calculation of the continued fraction for the I v I s, the Wronskian and the application of forward recurrence relations for the K v I s and backward recurrence for the I v I s. The main feature of the algorithm is that it does not require recalculations using normalization relations nor trial values to start the recurrences; the code evaluates in each step (already normalized) Bessel functions. The accuracy of the method (10 −16 for half-integral order and better than 2×10 −7 for integral order in our code) is limited only by the precision in the initial values for the recurrence and the maximum order available for a given value of the argument is restricted only by the maximum real number available in the computer.