  We use the Illustris-1 simulation to explore the capabilities of the
$\textit{Hubble}$ and $\textit{James Webb Space Telescope}$ data to analyze the
stellar populations in high-redshift galaxies, taking advantage of the combined
depth, spatial resolution, and wavelength coverage. For that purpose, we use
simulated broad-band ACS, WFC3 and NIRCam data and 2-dimensional stellar
population synthesis (2D-SPS) to derive the integrated star formation history
(SFH) of massive (M$_{\ast}>10^{10}\,$M$_{\odot}$) simulated galaxies at
$1<z<4$ that evolve into a local M$_{\ast}>10^{11}\,$M$_{\odot}$ galaxy. In
particular, we explore the potential of HST and JWST datasets reaching a depth
similar to those of the CANDELS and ongoing CEERS observations, respectively,
and concentrate on determining the capabilities of this dataset for
characterizing the first episodes in the SFH of local
M$_{\ast}>10^{11}\,$M$_{\odot}$ galaxies by studying their progenitors at
$z>1$. The 2D-SPS method presented in this paper has been calibrated to
robustly recover the cosmic times when the first star formation episodes
occurred in massive galaxies, i.e., the first stages in their integrated SFHs.
In particular, we discuss the times when the first 1% to 50% of their total
stellar mass formed in the simulation. We demonstrate that we can recover these
ages with typical median systematic offset of less than 5% and scatter around
20%-30%. According to our measurements on Illustris data, we are able to
recover that local M$_{\ast}>10^{11}\,$M$_{\odot}$ galaxies would have started
their formation by $z=16$, forming the first 5% of their stellar mass present
at $z \sim 1$ by $z=4.5$, 10% by $z=3.7$, and 25% by $z=2.7$.
