We present a new formulation to compute numerically stationary and axisymmetric equilibria of magnetized and self-gravitating astrophysical fluids. Under the assumption of ideal MHD, the stream function for the flow can be chosen as a basic variable with which the Euler-Maxwell equations are cast into a set of basic equations, i.e. a generalized Bernoulli equation and a Grad-Shafranov-like equation by employing various integral conditions. A novel feature of this formulation is that systems with stars, disks and winds are treated in a simple unified picture and the magnetic field structures can contain both poloidal and toroidal components.