  We present the first numerical solutions of the coupled Einstein-Maxwell
equations describing rapidly rotating neutron stars endowed with a magnetic
field. These solutions are fully relativistic and self-consistent, all the
effects of the electromagnetic field on the star's equilibrium (Lorentz force,
spacetime curvature generated by the electromagnetic stress-energy) being taken
into account. The magnetic field is axisymmetric and poloidal. Five dense
matter equations of state are employed. The partial differential equation
system is integrated by means of a pseudo-spectral method. Various tests passed
by the numerical code are presented. The effects of the magnetic field on
neutron stars structure are then investigated, especially by comparing
magnetized and non-magnetized configurations with the same baryon number. The
deformation of the star induced by the magnetic field is important only for
huge values of B (B>10^{10} T). The maximum mass as well as the maximum
rotational velocity are found to increase with the magnetic field. The maximum
allowable poloidal magnetic field is of the order of 10^{14} T (10^{18} G) and
is reached when the magnetic pressure is comparable to the fluid pressure at
the centre of the star. For such values, the maximum mass of neutron stars is
found to increase by 13 to 29 % (depending upon the EOS) with respect to the
maximum mass of non-magnetized stars.
