Metti il sistema di N equazioni in forma implicita:
f_1(x,y,z,...) = 0
...
f_N(x,y,z,..) = 0
poi definisci una funzione "fun" che ha come input il vettore delle tue variabili, e come output i valori di f_1,...,f_N calcolati su quelle variabili. A questo punto usi:
x = fsolve(fun,x0)
dove x0 è un vettore (x,y,z,...) attorno a cui vuoi cercare uno zero, cioè un vettore che risolve il sistema di equazioni. Però non so se usa Newton-Raphson.