We study solutions of the Grad-Shafranov equation which describe ideal magnetohydrodynamic (MHD) flows around a magnetized, rotating neutron star, where the rotation and magnetic axes are identical. The "force-free" limit of this equation is identical to the so-called pulsar equation. We derive the pulsar equation from a variational principle that minimizes the electromagnetic field energy subject to the constraints of fixed total angular momentum and total magnetic helicity. This procedure allows one to apply a result from fusion plasma theory, which suggests a specific form of the magnetic helicity for the configuration. This removes an indeterminancy from the pulsar equation which was not resolved in earlier attempts to find a unique equation for the electromagnetic fields. The pulsar equation is found to be a piecewise linear partial differential equation with a class of solutions that include self-collimated electromagnetic jets along the rotation axis of the star. These jets are confined by magnetic pinching and carry energy, angular momentum, and electric current far away from the star. We calculate numerically the global field structure for pulsar magnetospheres that contain electromagnetic jets. A principal result is that all of the poloidal magnetic field is confined within the light cylinder of the star, either as closed field loops connected to the surface of the star or as part of the jets. This yields a unique magnetospheric field topology with specific jet characteristics. We discuss the applicability of this model to compact synchrotron nebulae that have been observed in supernova remnants.