  We construct time-evolving gravitational potential models for a Milky
Way-mass galaxy from the FIRE-2 suite of cosmological-baryonic simulations
using basis function expansions. These models capture the angular variation
with spherical harmonics for the halo and azimuthal harmonics for the disk, and
the radial or meridional plane variation with splines. We fit low-order
expansions (4 angular/harmonic terms) to the galaxy's potential for each
snapshot, spaced roughly 25 Myr apart, over the last 4 Gyr of its evolution,
then extract the forces at discrete times and interpolate them between adjacent
snapshots for forward orbit integration. Our method reconstructs the forces
felt by simulation particles with high fidelity, with 95% of both stars and
dark matter, outside of self-gravitating subhalos, exhibiting errors $\leq4\%$
in both the disk and the halo. Imposing symmetry on the model systematically
increases these errors, particularly for disk particles, which show greater
sensitivity to imposed symmetries. The majority of orbits recovered using the
models exhibit positional errors $\leq10\%$ for 2-3 orbital periods, with
higher errors for orbits that spend more time near the galactic center.
Approximate integrals of motion are retrieved with high accuracy even with a
larger potential sampling interval of 200 Myr. After 4 Gyr of integration, 43%
and 70% of orbits have total energy and angular momentum errors within 10%,
respectively. Consequently, there is higher reliability in orbital shape
parameters such as pericenters and apocenters, with errors $\sim10\%$ even
after multiple orbital periods. These techniques have diverse applications,
including studying satellite disruption in cosmological contexts.
