I’m currently using Mathematica to symbolically solve a system of non-linear ODEs (10 equations with 8 variables) as I believe this is not yet possible using Symbolics.jl. Is anyone able to comment on whether it would be possible to call Mathematica into Julia using MathLink.jl to solve the below system. The documentation provided a simple example that I wasn’t able to extrapolate to my Mathematica code (below) successfully. Thanks very much.

The system is governed by 14 parameters:

```
kb1 = 0.01;
K1 = 10.0;
d1 = 1.0;
l11 = 1.0;
m11 = 0.01;
l12 = 1.0;
m12 = 5.0;
kb2 = 0.01;
K2 = 30.0;
d2 = 1.0;
l22 = 1.0;
m22 = 0.01;
l21 = 1.0;
m21 = 5.0;
```

and is solved in the following way in Mathematica:

```
Solve[{m11*X11 + m12*X12 - l11*M1*X10 - l12*M2*X10 == 0,
K1*X11 + kb1*X10 + m11*X11 + m21*X21 - d1*M1 - l11*M1*X10 - l21*M1*X20 == 0,
l11*M1*X10 - m11*X11 == 0,
m12*X12 + kb2*X20 + K2*X22 + m22*X22 - d2*M2 - l12*M2*X10 - l22*M2*X20 == 0,
l12*M2*X10 - m12*X12 == 0,
m21*X21 + m22*X22 - l21*M1*X20 - l22*M2*X20 == 0,
l22*M2*X20 - m22*X22 == 0, l21*M1*X20 - m21*X21 == 0,
X10 + X12 + X11 == 1,
X20 + X21 + X22 == 1
&& M1 >= 0
&& M2 >= 0} ,
{X10, M1, X11, M2, X12, X20, X22, X21}, Reals]
```

generating the below fixed points for variables M1 and M2:

```
{{X10 -> 0.000999999, M1 -> 9.93007, X11 -> 0.993006, M2 -> 29.9701,
X12 -> 0.00599402, X20 -> 0.000333333, X22 -> 0.999005,
X21 -> 0.000662005}}
```

How could I obtain this result in Julia?