Many astrophysical objects, including neutron stars and accretion disks, are powered by a combination of very strong magnetic fields and rotation.  The fields are so strong that magnetospheres of these objects can create electron-positron pairs, accelerate particles and result in broad range of radio to gamma-ray emission observed with high-energy missions and telescopes.  Our ability to understand these extreme environments has been hampered by the difficulty of solving the self-consistent behavior of strongly magnetized relativistic plasmas. I will describe recent progress in numerical modeling of magnetically-dominated plasmas and show applications to several sources of interest. I will present the numerical solution of the structure of pulsar magnetospheres, which has been a challenging problem for close to 40 years.   The knowledge of the magnetospheric shape together with the new observations of gamma-ray emission from pulsars with Fermi Telescope allow to directly constrain the location and physics of the acceleration regions in the magnetosphere. I will also discuss applications to magnetospheres of flaring magnetars and coronae of accretion disks, and conclude with comparison to observations and future directions.