I present our work on simulating the Glasma in the early stages of heavy ion collisions in a non-boost-invariant setting. Our simulation is based on the colored particle-in-cell method, which is used to numerically solve the Yang-Mills equations in 3+1 dimensions. This approach allows us to describe colliding nuclei with finite longitudinal width by extending the McLerran-Venugopalan model to include a parameter for the Lorentz-contracted but finite extent of the nucleus in the beam direction. We determine the rapidity profile of the Glasma energy density, which shows strong deviations from the boost-invariant result. Varying the parameters both broad and narrow profiles can be produced. We find reasonable agreement when we compare the results to rapidity profiles of measured pion multiplicities from RHIC.