We describe an algorithm to evaluate a wide class of functions and their derivatives, to extreme precision (25–30S) if required, which does not use any function calls other than square root. The functions are the Coulomb functions of positive argument ( F λ( x, η), G λ( x, η), x > 0, η, λ real) and hence, as special cases with η = 0, the cylindrical Bessel functions ( J μ( x), Y μ( x), x > 0, μ real), the spherical Bessel functions ( i λ( x), y λ( x), x > 0, λ real), Airy functions of negative argument Ai(- x), Bi(- x) and others. The present method has a number of attractive features: both the regular and irregular solution are calculated, all others of the functions can be produced from a specified minimum (not necessarily zero) to a specified maximum, functions of a single order can be found without all of the orders from zero, the derivatives of the functions arise naturally in the solution and are readily available, the results are available to different precisions from the same subroutine (in contrast to rational approximation techniques) and the methods can be used for estimating final accuracies. In addition, the sole constant required in the algorithm is π, no precalculated arrays of coefficients are needed, and the final accuracy is not dependent on that of other subroutines. The method works most efficiently in the region x ≈ 0.5 to x ≈ 1000 but outside this region the results are still reliable, even though the number of iterations within the subroutine rises. Even in these more asymptotic regions the unchanged algorithm can be used with known accuracy to test other specific subroutines more appropriate to these regions. The algorithm uses the recursion relations satisfied by the Coulomb functions and contains a significant advance over Miller's method for evaluating the ratio of successive minimal solutions ( F λ+1 / F λ ). It relies on the evaluation of two continued fractions and no infinite series is required for normalisation: instead the Wronskian is used.