Generalizing the inputs of the nlsolve function in Julia

This question has already been asked on another platform, but I haven't got an answer yet.

https://discourse.julialang.org/t/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)

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.

Two factors involved in a solution here. One is that NonlinearSolve.jl has explicit support for parameters, so you can just pass in parameters for NN and S and all. But now if you want to build equations symbolically, NonlinearSolve.jl comes with connections to ModelingToolkit.jl for generating numerical problems in an easy way. Here's an example of doing so:

using ModelingToolkit, NonlinearSolve

@variables x y z
@parameters σ ρ β

# Define a nonlinear system
eqs = [0 ~ σ * (y - x),
    0 ~ x * (ρ - z) - y,
    0 ~ x * y - β * z]
@named ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])

guess = [x => 1.0,
    y => 0.0,
    z => 0.0]

ps = [σ => 10.0
    ρ => 26.0
    β => 8 / 3]

prob = NonlinearProblem(ns, guess, ps)
sol = solve(prob, NewtonRaphson())

Thus to make use of this, you can simply change from SymPy to Symbolics expressions and use the MTK NonlinearSystem to construct the numerical problem from the symbolic expressions. On your problem this looks like:

using ModelingToolkit, NonlinearSolve
NN = 1
S = 1
@variables p[1:NN,1:S]
i = 1
j = 1
-331.330122303069*p[i,j]^(1.0) + p[i,j]^(2.81818181818182) -  1895.10478893046/(p[i,j]^(-1.0))^(2.0)

eqs = (@. 0 ~ -331.330122303069*p^(1.0) + p^(2.81818181818182) -  1895.10478893046/(p^(-1.0))^(2.0)) |> collect |> vec
@named sys = NonlinearSystem(eqs, vec(p), [])
prob = NonlinearProblem(sys, [p[1,1] => 1.0])
solve(prob)

#=

retcode: Success
u: 1-element Vector{Float64}:
 9.947428729753849e-23
 
=#

Note that we do not recommend NLsolve and have deprecated it for NonlinearSolve.jl throughout all of the numerical tools of SciML.