This paper studies the automation of the distortion analysis of multistage operational amplifiers by a symbolic analysis. The goal is to generate the readable closed-form design equations that relate the frequency-dependent harmonic distortion (HD) features, such as dc level, corner frequency, peaking, and zero, to the circuit parameters (including the nonlinearity coefficients). Prior research achieved such kind of characterization by carrying out the traditional signal-flow graph (SFG) analysis, which is highly algebra-involved and depends mainly on the manual derivation. This paper achieves automation by formulating a multiport circuit analysis strategy and applying a topological symbolic analysis method as an intermediate step in generating nodal transfer functions. The frequency-dependent HD factors are then characterized by analytically forming proper virtual current sources associated with the nonlinearity of each stage, modeled by a cubic polynomial. Implementation shows that the proposed tool can generate approximate closed-form HD characterization with comparable accuracy to full transistor-level simulation. The proposed automation tool can substantially mitigate heavy manual derivation work required by the existing SFG-based methods. Theory, software implementation, and circuit application are covered in detail.