We present a comprehensive analytic model of a relativistic jet propagation in expanding media. This model is the first to cover the entire jet evolution from early to late times, as well as a range of configurations that are relevant to binary neutron star mergers. These include low and high luminosity jets, unmagnetized and mildly magnetized jets, time-dependent luminosity jets, and Newtonian and relativistic head velocities. We also extend the existing solution of jets in a static medium to power-law density media with index $alpha<5$. Our model, which is tested and calibrated by a suite of 3D RMHD simulations, provides simple analytic formulae for the jet head propagation and breakout times, as well as a simple breakout criterion which depends only on the jet to ejecta energy ratio and jet opening angle. Assuming a delay time $ t_d $ between the onset of a homologous ejecta expansion and jet launching, the system evolution has two main regimes: strong and weak jets. The regime depends on the ratio between the jet head velocity in the ejecta frame and the local ejecta velocity, denoted as $ eta $. Strong jets start their propagation in the ejecta on a timescale shorter than $t_d$ with $eta gg 1$, and within several ejecta dynamical times $eta$ drops below unity. Weak jets are unable to penetrate the ejecta at first (start with $eta ll 1$), and breach the ejecta only after the ejecta expands over a timescale longer than $ t_d $, thus their evolution is independent of $ t_d $. After enough time, both strong and weak jets approach an asymptotic phase where $eta$ is constant. Applying our model to short GRBs, we find that there is most likely a large diversity of ejecta mass, where mass $ lesssim 10^{-3}~{rm M}_{odot} $ (at least along the poles) is common.