Generalizing the inputs of the nlsolve function in Julia

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)


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

Now, I find the zeros

function K(F,p)
F[1] = resf2(p[1])

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.