We construct equilibrium configurations of magnetized, two-fluid neutron stars using an iterative numerical method. Working in Newtonian framework we assume that the neutron star has two regions: the core, which is modelled as a two-component fluid consisting of type-II superconducting protons and superfluid neutrons, and the crust, a region composed of normal matter. Taking a new step towards more complete equilibrium models, we include the effect of entrainment, which implies that a magnetic force acts on neutrons, too. We consider purely poloidal field cases and present improvements to an earlier numerical scheme for solving equilibrium equations, by introducing new convergence criteria. We find that entrainment results in qualitative differences in the structure of field lines along the magnetic axis.