I have two functions $P_1,Q_1$ which are both depending on the variables $r_{T_m}, r_{i_f}$, which are theirselves functions of $T_m, mi_f$. I am interested in computing the Jacobian

$$frac{partial (P_1,Q_1)}{partial(r_{T_m}^2,r_{i_f}^2)}Bigg|_{r_{Tm}=f(T_m) \ r_{i_f}=g(mi_f)}.$$

To do that, I first define the function

```
powerD(f_, x_^(k_.)) := powerD(f, {x^k, 1});
powerD(f_, {x_^(k_.), 0}) := f;
powerD(f_, vars__) := Fold(powerD, f, {vars});
powerD(f_, {x_^(k_.), n_Integer?Positive}) := Det(Append(Table((j!/i!) Binomial(k i, j) x^(k i - j), {i, n - 1}, {j, n}), Table(D(f, {x, j}), {j, n})))/(k x^(k - 1))^Binomial(n + 1, 2);
```

(because I want to compute the derivative wrt the power of $r_{T_m}, r_{mi_f}$), then I want to substitute $r_{T_m}=f(T_m)$ and $r_{i_f}=g(mi_f)$. However, something in the code fails at this point. I think it is because instead of first computing the partial derivatives $frac{partial P_1}{partial r^2_{T_m}}$ etc. and then substituting $r_{T_m}=f(T_m)$, it tries to compute $frac{partial P_1}{partial f^2(T_m)}$.

I attach here the code ($P_1,Q_1$ are actually functions of other parameters too, but I substitute these other parameters with the needed numerical values before computing the Jacobian).

The function $P_1(r_{T_m},mi_f)$ is given by

```
P1(230 1.73, 0.2, 0.0072, 100 3.14, rTm, rif)=-1.12469*10^-8 (1.71055*10^13 - 112.657 rif^2 + 112.657 rTm^2 - 9.60432*10^-10 (-5.66361*10^20 - 3.52759*10^9 rif^2 + 3.52759*10^9 rTm^2 + Sqrt((5.66361*10^20 + 3.52759*10^9 rif^2 - 3.52759*10^9 rTm^2)^2 - 3.26026*10^16 (3.2715*10^26 - 3.83021*10^15 rif^2 + 12691.7 rif^4 - 4.07533*10^15 rTm^2 - 25383.3 rif^2 rTm^2 + 12691.7 rTm^4))))
```

and the function $Q_1(r_{T_m},mi_f)$ is given by

```
Q1(230 1.73, 0.2, 0.0072, 100 3.14, rTm, rif)=6.13447*10^-17 (-5.66361*10^20 - 3.52759*10^9 rif^2 + 3.52759*10^9 rTm^2 + Sqrt((5.66361*10^20 + 3.52759*10^9 rif^2 - 3.52759*10^9 rTm^2)^2 - 3.26026*10^16 (3.2715*10^26 - 3.83021*10^15 rif^2 + 12691.7 rif^4 - 4.07533*10^15 rTm^2 - 25383.3 rif^2 rTm^2 + 12691.7 rTm^4)))
```

I compute the elements of the (2×2) Jacobian matrix as

```
J11(rTm_, rif_) := FullSimplify(powerD(P1(230 1.73, 0.2, 0.0072, 100 3.14, rTm, rif), rTm^2))
J12(rTm_, rif_) := FullSimplify(powerD(P1(230 1.73, 0.2, 0.0072, 100 3.14, rTm, rif), rif^2))
J21(rTm_, rif_) := FullSimplify(powerD(Q1(230 1.73, 0.2, 0.0072, 100 3.14, rTm, rif), rTm^2))
J22(rTm_, rif_) := FullSimplify(powerD(Q1(230 1.73, 0.2, 0.0072, 100 3.14, rTm, rif), rif^2))
```

and then I build the full Jacobian matrix as

```
J(rTm_, rif_) := {{J11(rTm, rif), J12(rTm, rif)}, {J21(rTm, rif),J22(rTm, rif)}}
```

Up to here everything is okay, but then when I substitute

$$r_{T_m}=f(T_m)=sqrtfrac{(230cdot1.73)^4 + 4cdot(230cdot1.73)^2cdot0.2cdot100cdot3.14cdot Tm}{4cdot(0.2)^2}$$

and

$$r_{i_f}=g(mi_f)=sqrtfrac{230cdot1.73cdot mifcdot100cdot3.14 }{(0.2)^2 + (100cdot3.14cdot 0.0072)^2}$$

with the command

```
J(Sqrt(((230 1.73)^4 + 4 (230 1.73)^2 0.2 100 3.14 Tm)/(4 (0.2)^2)), Sqrt(((230 1.73) mif 100 3.14)/((0.2)^2 + (100 3.14 0.0072)^2)))
```

I get the errors

“6.25 (2.50666*10^10+3.97711*10^7 Tm) is not a valid variable”,

“24254.6mif is not a valid variable.”

I have just started today using the software, so I am sure it is a pretty easy mistake to be solved, maybe I should assign the values differently to the various Jacobian elements $J_{i,j}$.