Here's the answer I found out:
The analytical equations can be solved as a convex optimization problem. i.e by keeping all terms on one side and introducing an inequality constraint. For example, in this problem,
dP/drho = (dRp/drho)*Rp' + Rp*(dRp/drho)'
can be written as
((dRp/drho)*Rp' + Rp*(dRp/drho)' -dP/drho) >=0
and then solve an optimiation ...