After an extensive process usyng the SymPy in Julia, I generated a system of nonlinear equations. My system is allocated in a matrix NxS. Something like this(NN = 2, S = 2).

I would like to adapt the system to use the NLsolve package. I do some boondoggle for the case NN=1 and S =1. The system_equations2 function give me the nonlinear system, like the figure

```
using SymPy
using Plots
using NLsolve
res = system_equations2()
```

In order to simulate the output, I do this:

```
NN = 1
S = 1
p= [Sym("p$i$j") for i in 1:NN,j in 1:S]
res = [ Eq( -331.330122303069*p[i,j]^(1.0) + p[i,j]^(2.81818181818182) - 1895.10478893046/(p[i,j]^(-1.0))^(2.0),0 ) for i in 1:NN,j in 1:S]
resf = convert( Function, lhs( res[1,1] ) )
plot(resf, 0 ,10731)
```

Now

```
resf = convert( Function, lhs( res[1,1] ) )
# This for the argument in the nlsolve function
function resf2(p)
p = Tuple(p)[1]
r = resf(p)
return r
end
```

Now, I find the zeros

```
function K(F,p)
F[1] = resf2(p[1])
end
nlsolve(K , [7500.8])
```

I would like to generalize this price to any NN and any S. I believe there is a simpler way to do this.