A Maple algorithm for the computation of the zeros of orthogonal polynomials (OPs) and special functions (SFs) in a given interval [x1, x2] is presented. The program combines symbolic and numerical calculations and it is based on fixed point iterations. The program uses as inputs the analytic expressions for the coefficients of the three-term recurrence relation and a differencedifferential relation satisfied by the set of OPs or SFs. The performance of the method is illustrated with several examples: Hermite, Chebyshev, Legendre, Jacobi and Gegenbauer polynomials, Bessel, Coulomb and Conical functions.