We present our results on genuinely 3+1D simulations of the Glasma in heavy ion collisions. By incorporating nuclei with finite thickness along the beam axis, we naturally break boost invariance and obtain Gaussian-like rapidity profiles of energy densities already at the tree level . The profiles resemble strong coupling results and agree surprisingly well with experimental data of pion multiplicities as obtained at RHIC. The applicability of these simulations is limited by the numerical Cherenkov instability. We construct a new lattice action which allows us to cure this instability along one lattice direction . This paves the way for simulations at higher energies and better resolution.