Using a new formulation to compute structures of stationary and axisymmetric magnetized barotropic stars in Newtonian gravity, we have succeeded in obtaining numerically exact models of stars with extremely high interior magnetic fields. In this formulation, there appear four arbitrary functions of the magnetic flux function from the integrability conditions among the basic equations. Since in our new formulation these arbitrary functions appear in the expression of the current density, configurations with different current distributions can be specified by choosing the forms of the arbitrary functions.

By choosing appropriate forms for the four arbitrary functions, we have solved many kinds of equilibrium configurations both with poloidal and toroidal magnetic fields. Among them, by choosing special form for the toroidal current density, we have been able to obtain magnetized stars which have extremely strong poloidal magnetic fields deep inside the core region near the symmetric axis. By adopting the appropriate model parameters for the neutron stars, the magnetic fields could be 1014 ~ 1015 G on the surfaces and be about 1017 G in the deep interior regions. For other model parameters appropriate for white dwarfs, the magnetic fields could be around 107 ~ 108 G (surface regions) and 109 ~ 1010 G (core regions). It is remarkable that the regions with very strong interior magnetic fields are confined to a very narrow region around the symmetric axis in the central part of the stars. The issues of stability of these configurations and of evolutionary paths to reach such configurations need to be investigated in the future work.