This work is a sequel to our two 2023 publications [PLB 837 137669,
NPA 1039 122725] where fitting 14 1$s_\Lambda$ and 1$p_\Lambda$
single-particle binding energies in hypernuclei across the periodic table led
to a well-defined $\Lambda$-nucleus optical potential. The potential consists
of a Pauli modified linear-density ($\Lambda N$) and a quadratic-density
($\Lambda NN$) terms. The present work reports on extending the above analysis
to 21 $\Lambda$ single-particle data points input by including 1$d_\Lambda$
and 1$f_\Lambda$ states in medium-weight and heavy hypernuclei.
The upgraded results for the $\Lambda N$ and $\Lambda NN$
potential depths at nuclear-matter density $\rho_0=0.17$~fm$^{-3}$,
$D^{(2)}_\Lambda=-37.5\mp 0.7$~MeV and $D^{(3)}_\Lambda=+9.8\pm 1.2$~MeV
together with the total depth $D_\Lambda=-27.7\pm 0.5$~MeV, agree within
errors with the earlier results. The $\Lambda$ hypernuclear overbinding
associated with the $\Lambda N$-induced potential depth $D^{(2)}_\Lambda$
agrees quantitatively with a recent combined analysis of low-energy
$\Lambda p$ scattering data and correlation functions [PLB 850 (2024) 138550].
These results, particularly the size of the repulsive $D^{(3)}_\Lambda$,
provide an essential input towards resolving the 'hyperon puzzle' in the core
of neutron stars. We also show that a key property of our $\Lambda NN$-induced
potential term, i.e. a need to suppress the quadratic-density $\Lambda NN$
term involving an excess neutron and a $N=Z$ core nucleon, can be tested in
the forthcoming JLab E12-15-008 experiment.