In order to compute the free core nutation of the terrestrial planets, such as Earth and Mars, the Moon and lower degree normal modes of the Jovian planets, we propose a linear operator method (LOM). Generalized surface spherical harmonics (GSSHs) are usually applied to the elliptical models with a stress tensor, which cannot be expressed in vector spherical harmonics explicitly. However, GSSHs involve complicated math. LOM is an alternative to GSSHs, whereas it only deals with the coupling fields of the same azimuthal order m, as is the case when a planet model is axially symmetric and rotates about that symmetry axis. We extend LOM to any asymmetric 3D model. The lower degree spheroidal modes of the Earth are computed to validate our method, and the results agree very well with what is observed. We also compute the normal modes of a two-layer Saturn model as a simple application.