A scheme for the computation of the zeros of special functions and orthogonal polynomials is developed. The structure of the first order difference-differential equations (DDEs) satisfied by two fundamental sets of solutions of second order ODEs
being the order of the solutions and
a family of continuous functions is studied. The structure of the ODEs is also used to set global bounds on the differences between adjacent zeros of functions of consecutive orders and to find iteration steps which guarantee that all the zeros inside a given interval can be found with certainty. As illustration is described how to implement this sequence to the calculation of the zeros of arbitrary solutions of certain known equations (Bessel, Coulomb, Legendre, Hermite, Laguerre).