We present results of the first self-consistent, time-dependent, 2-D calculations of the accretion of plasma onto polar caps of high luminosity (L*>1036erg-s−1) magnetized neutron stars. We follow the temporal and spatial evolution of three fluids, electrons, ions and photons in a superstrong (B=3×1012 Gauss) dipole magnetic field where radiation pressure dominates plasma pressure by solving coupled 2-D equations of radiation hydrodynamics. We have included several physical processes in the radiation-plasma coupling in superstrong magnetic fields (Klein, et al., 1984, Santa Cruz Workshop on High Energy Transients, and Arons, this conference). We solve the resulting system of coupled 2-D PDEs on a Cray XMP-48 by applying implicit finite-difference techniques with iterative operator splitting methods. We present results for two models of 5×1037 erg-s and 1.5×1038 erg-s−1 super-Eddington luminosity on one polar cap, each having initial mass flux independent of co-latitude of a field lines footprint. We find (a) Radiation develops a broad transverse fan beam that emerges from an annulus 0.2–0.5km above the polar cap. (b) The beam profile is determined by advective trapping of radiation in optically thick (τ11,τ⊥ ≈103) flow. Here the time for diffusion of radiation up through the accretion column is ≫ the time for downward advection. (c) There is a three fluid nonequilibrium with Ti≫Tγ≥Te. (d) Maximum photon temperature of ≈ 10–20 keV in the fan beam is in the observed range. (e) Cyclotron emission ≫ bremsstrahlung as a source of photons. (f) At early times (≪lms) radiation pressure strongly decelerates flow to 10−3 of freefall in central regions of accretion column resulting in a density mound, but plasma freefalls down the sides of the column. (g) Analytical models have reasonable agreement with numerical calculations; velocity and energy density roughly Gaussian transversally and exponential vertically, until the onset of “photon bubbles” after several dynamical times (∼lms). (h) Multiple “photon bubbles” rising subsonically in the accretion column form in the high luminosity model. We believe the photon bubbles to be a possible consequence of overstable convection in super-Eddington flows. These photon bubbles could be observable as 10–100μs fluctuations in the emergent flux and, thus, be an important diagnostic for inhomogeneous structure of the column.