The singularity structure of a second-order ordinary differential equation with polynomial coefficients often yields the type of solution. If the solution is a special function that is studied in the literature, then the result is more manageable using the properties of that function. It is straightforward to find the regular and irregular singular points of such an equation by a computer algebra system. However, one needs the corresponding indices for a full analysis of the singularity structure. It is shown that the $theta$-operator method can be used as a symbolic computational approach to obtain the indicial equation and the recurrence relation. Consequently, the singularity structure which can be visualized through a Riemann P-symbol leads to the transformations that yield a solution in terms of a special function, if the equation is suitable. Hypergeometric and Heun-type equations are mostly employed in physical applications. Thus only these equations and their confluent types are considered with SageMath routines which are assembled in the open-source package symODE2.