I didn’t find a way to deal with vector using pkg NLSolve in the following problem.

Let q=[x,y] and p=[u,v] and defining the functions f and g as f(q) = q / (√(q'*q))^3,\ g(p) = p. Now, giving q_0=[0.4,0] and p_0=[0,0.5]. I need to solve by Newton-Raphson method the following

Your parameters can’t be a vector of vectors, you need to do something like

ini_gauss = [q₀..., p₀...]

i.e. you need to flatten this structure. You’ll also need to adjust that closure, i.e.

(q,p)-> F(q₀,p₀,q,p)

should be

qp -> F(q₀,p₀, view(qp, 1:2), view(qp, 3:4))

and similarly, F should return a flat vector instead of a vector of vectors:

function F(q₀,p₀,q,p)
p′ = q - p₀ + 0.5 * f((q₀+q)/2)
q′ = p - q₀ - 0.5 * g((p₀+p)/2)
[p′..., q′...]
end

Note, the above code is not very optimized and there’s a lot that could be done to make it faster. However, at least with the initial conditions you gave, it seems this system of equations does not converge to a solution with any of the solvers I tried, so I wonder if there’s a problem with your initial conditions or your equations.