This work presents an analytical method to derive the stress fields for an elastic medium containing a circular anisotropic inclusion embedded in an infinitely extended isotropic matrix subjected to far-field uniform stresses using symbolic software. The Airy stress functions are employed to represent the stress distribution in such materials. Traditionally, the Airy stress function has been mainly used to derive the stress fields for 2-D isotropic materials. For 2-D anisotropic materials, Lekhnitskii’s formulation provides an alternative approach, where the stress function is expressed as the real part of the sum of two analytic functions of two independent complex variables, determined from the roots of the characteristic equation associated with the anisotropic elastic constants. In this study, the stress function for isotropic materials is used for the matrix phase, while the stress function by Lekhnitskii is used for the inclusion phase. The stress function for the matrix phase is expressed by the Laurent series, whereas the Taylor series is employed for the inclusion phase. Determining the unknown complex coefficients of these series requires satisfying the continuity conditions of the displacement and traction force across the interface between the matrix and inclusion. The obtained results are novel and in closed form. The proposed method can be extended to other problems, such as thermal stress and heat conduction, demonstrating its broad applicability across various fields.