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 α < 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 η. Strong jets start their propagation in the ejecta on a timescale shorter than t d with η 1, and within several ejecta dynamical times η drops below unity. Weak jets are unable to penetrate the ejecta at first (start with η 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 η is constant. Applying our model to short GRBs, we find that there is most likely a large diversity of ejecta mass, where mass 10 −3 M (at least along the poles) is common.