We present population-orbit superposition models for external galaxies based on Schwarzschild's orbit-superposition method, by tagging the orbits with age and metallicity. The models fit the density distributions, kinematic, and age and metallicity maps from integral field unit (IFU) spectroscopy observations. We validate the method and demonstrate its power by applying it to mock data, similar to those obtained by the Multi-Unit Spectroscopic Explorer (MUSE) IFU on the Very Large Telescope (VLT). These mock data are created from Auriga galaxy simulations, viewed at three different inclination angles (θ = 40°, 60°, 80°). Constrained by MUSE-like mock data, our model can recover the galaxy's stellar orbit distribution projected in orbital circularity λz versus radius r, the intrinsic stellar population distribution in age t versus metallicity Z, and the correlation between orbits' circularity λz and stellar age t. A physically motivated age─metallicity relation improves the recovering of intrinsic stellar population distributions. We decompose galaxies into cold, warm, and hot+counter-rotating components based on their orbit circularity distribution, and find that the surface density, velocity, velocity dispersion, and age and metallicity maps of each component from our models well reproduce those from simulation, especially for projections close to edge-on. These galaxies exhibit strong global age versus σz relation, which is well recovered by our model. The method has the power to reveal the detailed build-up of stellar structures in galaxies, and offers a complement to local resolved, and high-redshift studies of galaxy evolution.