  In this work, we analyse the density profiles of subhaloes with masses
$M_\mathrm{sh} \geq 1.4 \times 10^8$ M$_\odot$ in the TNG50 simulation, with
the aim of including baryonic effects. We evaluate the performance of
frequently used models, such as the standard NFW, the Einasto, and a smoothly
truncated version of the NFW profile. We find that these models do not perform
well for the majority of subhaloes, with the NFW profile giving the worst fit
in most cases. This is primarily due to mismatches in the inner and outer
logarithmic slopes, which are significantly steeper for a large number of
subhaloes in the presence of baryons. To address this issue, we propose new
three-parameter models and show that they significantly improve the goodness of
fit independently of the subhalo's specific properties. Our best-performing
model is a modified version of the NFW profile with an inner log-slope of -2
and a variable truncation that is sharper and steeper than the slope transition
in the standard NFW profile. Additionally, we investigate how both the
parameter values of the best density profile model and the average density
profiles vary with subhalo mass, $V_\mathrm{max}$, distance from the host halo
centre, baryon content and infall time, and we also present explicit scaling
relations for the mean parameters of the individual profiles. The newly
proposed fit and the scaling relations are useful to predict the properties of
realistic subhaloes in the mass range $10^8$ M$_\odot$ $\leq M_\mathrm{sh}
\leq$ $10^{13}$ M$_\odot$ that can be influenced by the presence of baryons.
