That is not going to work in general. You run the risk of accidentally substituting into something that happens to match the value of diff(p(rho,m1,E),rho) but which has a different origin.
For example look at the term (E+p)*m1/rho and imagine it differentiated with regards to m1. That is going to give you E/rho + p(rho,m1,E)*rho + diff(p(rho,m1,E),m1)*m1/rho if I mentally constructed the chain rule correctly in the early hour of the morning. Now you take diff(p(rho,m1,E),rho) and want to substitute DpDrho for that. But suppose that happens to come out as E/rho ? Then you would substitute gamma*T for the E/rho even though the sub-expression was not at all generated by differentiation of p with respect to rho.
You are trying to "tag" the results of the partial differentiations to show their source, wanting to only do the substitution of DpDrho in places that arose from differentiating p with respect to rho. But the result of the differentiation is a value and the value might happen to have arisen independently through other means.
You also need to consider the possibility that the differentiation of p with respect to rho has an additive constant such as 59 and that when the overall expression that p appears in in F was differentiated with respect to rho, that the additive constant cancels out or modifies. X+59 might now show up as X+63, but still be sourced from the partial derivative and the X+63 should be substituted as (DpDrho +4). But subs() does pattern matching, not algebraic substitution, and if the pattern does not occur exactly as-is then the substitution will fail.
This method is just too weak and too undefined. You want a procedure that would not change if p was known before the jacobian is taken, but you are going to have difficulty constructing such a method as accidental linkages can change what the jacobian should come out as, under a series of conditions.
