Calculating Love numbers (tidal deformation) for non-homogeneous planets
I actually just found an article that explains such a calculation in detail, using a matrix propagator method: Padovan et al. 2018
This is a calculation for gas giants or bodies in hydrostatic equilibrium. The calculation of the Love numbers for terrestrial planets that have viscoelastic behaviour requires the rheology (which includes models for the rigidity). I did not find the equations that need to be solved in that case, but a similar matrix propagator method can be used apparently.
I will not repeat the matrix propagation method here, as it is quite well elaborated in the paper I just mentioned, except for the equation you have to start from, which is explained in another article, but looks a bit too complicated for me at this point.