Neutron stars are among the most fascinating astrophysical sources, being characterized by strong gravity, densities about the nuclear one or even above, and huge magnetic fields. Their observational signatures can be extremely diverse across the electromagnetic spectrum, ranging from the periodic and low-frequency signals of radio pulsars, up to the abrupt high-energy gamma-ray flares of magnetars, where energies of are released in a few seconds. Fast-rotating and highly magnetized neutron stars are expected to launch powerful relativistic winds, whose interaction with the supernova remnants gives rise to the non-thermal emission of pulsar wind nebulae, which are known cosmic accelerators of electrons and positrons up to PeV energies. In the extreme cases of proto-magnetars (magnetic fields of G and millisecond periods), a similar mechanism is likely to provide a viable engine for the still mysterious gamma-ray bursts. The key ingredient in all these spectacular manifestations of neutron stars is the presence of strong magnetic fields in their constituent plasma. Here we will present recent updates of a couple of state-of-the-art numerical investigations by the high-energy astrophysics group in Arcetri: a comprehensive modeling of the steady-state axisymmetric structure of rotating magnetized neutron stars in general relativity, and dynamical 3D MHD simulations of relativistic pulsar winds and their associated nebulae.