In this work, we develop an analytic treatment that would allow us to model the magnetic deformation of neutron stars in the quasi-Newtonian framework. We analyze the problem with taking density profiles from various analytic approaches. The first density profile is the n=1 solution of the Lane-Emden equation, the second one is a parabolic ansatz for density and the last density distribution is a fourth-order polynomial fit to a relativistic solution we obtained with a microphysical EoS. The last density profile is an accurate modeling obtained by fitting to the true density profile acquired numerically for various EoSs  and it allows us to have an analytical solution for the magnetic field profile for almost a universal EoS. Contour plots of purely poloidal field indicate that the poloidal ﬁeld components contribute most of the interiors of the star and have the maximum and minimum values of ﬁeld strength at the origin and the edge of the star, respectively. While the projection of the toroidal magnetic ﬁeld lines on the meridional plane shows that the toroidal ﬁeld distribute in wide regions in the vicinity of the equatorial plane and the maximum toroidal fields occur in the center of torus. We almost incorporate all the magnetic field configurations assessed from different density profiles. The excess mass calculated in this quasi-Newtonian framework is analogous to the excess mass obtained in the general relativity framework and the polynomial magnetic field profile [2,3].
 N. Jiang and K. Yagi, Phys. Rev. D 99, 124029 (2019).
 M. Zamani, and M. Bigdeli, J. Phys. G: Nucl. Part. Phys. 46, 075201 (2019)
 M. Zamani, and M. Bigdeli, Astronomische Nachrichten 342, 633 (2021)