We present a unified approach to calculating the zeros of the classical orthogonal polynomials and discuss the electrostatic interpretation and its connection to the energy minimization problem. This approach works for the generalized Bessel polynomials, including the normalized reversed variant, as well as the Viet\'e--Pell and Viet\'e--Pell--Lucas polynomials. We briefly discuss the electrostatic interpretation for each aforesaid case and some recent advances. We provide zeros and error estimates for various cases of the Jacobi, Hermite, and Laguerre polynomials and offer a brief discussion of how the method was implemented symbolically and numerically with Maple. In conclusion, we provide possible avenues for future research.