Note that your problem has precisely the form z = g(z) if you were to consider z := (x, y) and g := (g1, g2).
The only thing you have to do is to concatenate both variables and functions. Assuming x and y have the same size and dimensions, one way to do this would be
z = cat(x, y; dims=ndims(x)+1) # concatenate along a new dimension
And then a similar procedure for your functions.
But personally I am not a fan of this because explicit concatenation with these cat functions is not very efficient, especially for big arrays like yours. It also only works properly if x and y have the same size and typeof. On top of that, you’d need to throw @views everywhere to avoid further allocations.
So something much nicer would be to use an array of arrays such as
z = [x, y]
g(z) = [g1(z[1], z[2]), g2(z[1], z[2])]
But nlsolve can’t really operate on such containers.
However, it can work on slightly more refined containers such as the ones in RecursiveArraysTools, which are designed specifically for these tasks.
using RecursiveArrayTools
z = VectorOfArray([x, y])
g(z) = VectorOfArray([g1(z[1], z[2]), g2(z[1], z[2])])
nlsolve(g, z; method=:anderson)