Quiet@Remove["Global`*"];
$HistoryLength=2;
SetDirectory[NotebookDirectory[]];
mf=MatrixForm;
randomMatrix[rnm_]:=Table[RandomReal[{-1,1}]&[i],{i,rnm},{j,rnm}];

(*http://kfe.fjfi.cvut.cz/~limpouch/numet/nr/node14.html*)
(*http://kfe.fjfi.cvut.cz/~limpouch/numet/nr/node15.html*)
(*http://fourier.eng.hmc.edu/e161/lectures/ica/node13.html*)

randomMatrix[3]//mf

Map[g,{x1,x2,x3}]

g/@{x1,x2,x3}

(*xx={x1,x2,x3};
res={};
For[i=1,i<=Length[xx],i++,res=Append[res,g[xx[[i]]]]];
res*)

Range[n]

SeedRandom[1]
n=3;
{a,b,c}=randomMatrix[n]&/@Range[n];
{p1,p2,p3}=-3Random[]&/@Range[n]

ClearAll[x]
f1[xvect_]:=Total[(a.xvect)]+p1;
f2[xvect_]:=Plus@@(b.xvect)+p2;
f3[xvect_]:=Plus@@(c.xvect)+p3;
vars=(Subscript[x, #])&/@Range[n];
vectfces=Through[{f1,f2,f3}[vars]];
origEqns=Thread[vectfces==(0&/@Range[n])];

vectfces//mf
origEqns//mf

(Subscript[x, #])&/@Range[n]

(Subscript[x, #])&["Monkey"]

vars

fce2[x]:=x^2

fce2[x]
fce2[2]

?fce

x=5
fce[x_]:=x^2

fce[2]

fce3=Function[#^2]

fce3[2]

fce/@{2,4}
(Function[#^2])/@{2,4}
(#^2&)/@{2,4}

Function[#^2]

a

grad[expr_]:=D[expr,#]&/@vars;
jacobian=grad/@vectfces;
xStart=Function[Random[]]/@Range[n]
xNext[xPrevious_]:=xPrevious-Inverse[jacobian/.Thread[vars->xPrevious]].(vectfces/.Thread[vars->xPrevious]);
xNext[xStart]
Solve[origEqns]
(*for linear equations one step is enough*)

jacobian

jacobian

vectfces

grad[-2.0060780335801036`+0.09487851477344345` (Subscript[x, 1])-1.8321298293501713` (Subscript[x, 2])-0.4974583286175607` (Subscript[x, 3])]

grad/@vectfces

vectfces

vars

Total[vars]

grad[(Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])]

vars

grad[(Subscript[x, 1])+2*(Subscript[x, 2])]

(*nonlinear example, quadratic forms zeros*)
SeedRandom[1]
n=3;
{a,b,c}=randomMatrix[n]&/@Range[n];
{p1,p2,p3}=-3Random[]&/@Range[n];
f1[xvect_]:=(a.xvect).xvect+p1;
f2[xvect_]:=(b.xvect).xvect+p2;
f3[xvect_]:=(c.xvect).xvect+p3;
vars=(Subscript[x, #])&/@Range[n];
vectfces=Through[{f1,f2,f3}[vars]];
origEqns=Thread[vectfces==(0&/@Range[n])];
grad[expr_]:=D[expr,#]&/@vars;
jacobian=grad/@vectfces;
xStart=Random[]&/@Range[n]
xNext[xPrevious_]:=xPrevious-Inverse[jacobian/.Thread[vars->xPrevious]].(vectfces/.Thread[vars->xPrevious]);
xNext[xStart]
NestList[xNext,xStart,50][[-1]]
(*FindRoot[origEqns,Thread[List[vars,xStart]]]
Solve[origEqns]
*)

FindRoot[origEqns,Thread[List[vars,xStart]]]
Solve[origEqns]

vectfces

(*load flow problem, cartesian form of complex numbers*)
Quiet@Remove["Global`*"];
simplify=#&;(*"lower" programing languages have no Mathematica Simplify[#]&*)
(*simplify=Simplify;*)
conjugate[{re_,im_}]:={re,-im};
times[ua:{uare_,uaim_},ub:{ubre_,ubim_}]:={uare*ubre-uaim*ubim,uare*ubim+ubre*uaim};
re[{re_,im_}]:={re,0};
im[{re_,im_}]:={im,0};
sqabs[{re_,im_}]:={(re^2)+(im^2),0};
powerMinus1[z:{re_,im_}]:=(({re,-im})/((re^2)+(im^2)));
ratio[a_,b_]:=times[powerMinus1[b],a]
apparentPower[u_,i_]:=times[u,conjugate[i]];
activePower[u_,i_]:=re[apparentPower[u,i]];
reactivePower[u_,i_]:=im[apparentPower[u,i]];
current[u1_,u2_,z_]:=ratio[u1-u2,z];
toMyForm[z_]:={Re[z],Im[z]};
fromMyForm[{re_,im_}]:=re+I*im;
rearrangeRealEqn[eqn_Equal]:={simplify[eqn]}/.Equal[a_,b_]:>Equal[a-b,0];
rearrangeComplexEqn[{rea_,ima_}=={reb_,imb_}]:=DeleteCases[#,rearrangeRealEqn[True]]&@Flatten[rearrangeRealEqn/@{simplify[rea==reb],simplify[ima==imb]}];
rearrangeComplexEqn[{a,b}=={c,d}];
rearrangeRealEqn[x+2y==y-x];

kz=2.5;(*respects length of lines*)
kReg=2.5;(*ratio of powers of sources*)
{(Subscript[z, 1]),(Subscript[z, 2]),(Subscript[z, 3])}=kz*(toMyForm/@{.01+0.05I,2*(01+0.05I),0.01+0.02I});
powerFactor4=0.8;
powerFactor5=0.7;
S4=(10^3)*{powerFactor4,Sin[ArcCos[powerFactor4]]};
S5=500*{powerFactor5,Sin[ArcCos[powerFactor5]]};
{(Subscript[u, #])={(Subscript[ure, #]),(Subscript[uim, #])},(Subscript[i, #])={(Subscript[ire, #]),(Subscript[iim, #])}}&/@Range[5];
Unom=230.;
eqnVoltage3=rearrangeComplexEqn[sqabs[(Subscript[u, 3])]==toMyForm[(Unom^2)]];
eqnZ1=rearrangeComplexEqn[(Subscript[i, 1])==current[(Subscript[u, 1]),(Subscript[u, 3]),(Subscript[z, 1])]];
eqnZ2=rearrangeComplexEqn[(Subscript[i, 2])==current[(Subscript[u, 2]),(Subscript[u, 3]),(Subscript[z, 2])]];
eqnZ3=rearrangeComplexEqn[(Subscript[i, 3])==current[(Subscript[u, 3]),(Subscript[u, 4]),(Subscript[z, 3])]];
eqnNod3=rearrangeComplexEqn[(Subscript[i, 1])+(Subscript[i, 2])==(Subscript[i, 3])];
eqnNod4=rearrangeComplexEqn[(Subscript[i, 4])+(Subscript[i, 5])==(Subscript[i, 3])];
eqnSourceActivePowerRatio=rearrangeComplexEqn[activePower[(Subscript[u, 1]),(Subscript[i, 1])]==kReg*activePower[(Subscript[u, 2]),(Subscript[i, 2])]];
eqnSourcReactivePowerRatio=rearrangeComplexEqn[reactivePower[(Subscript[u, 1]),(Subscript[i, 1])]==kReg*reactivePower[(Subscript[u, 2]),(Subscript[i, 2])]];
eqnApparentPower4=rearrangeComplexEqn[apparentPower[(Subscript[u, 4]),(Subscript[i, 4])]==S4];
eqnApparentPower5=rearrangeComplexEqn[apparentPower[(Subscript[u, 4]),(Subscript[i, 5])]==S5];
eqnSettingAbsoluteTime={(Subscript[ure, 3])==0};
eqnsTogether=Union[eqnVoltage3,eqnZ1,eqnZ2,eqnZ3,eqnNod3,eqnNod4,eqnSourceActivePowerRatio,eqnSourcReactivePowerRatio,eqnApparentPower4,eqnApparentPower5,eqnSettingAbsoluteTime];
unknownVariables=Union@Cases[eqnsTogether,Subscript[_,_],{0,∞}];
Length/@{eqnsTogether,unknownVariables}

secondParamteterOfFindRoot=unknownVariables/.{(Subscript[iim, index_]):>{(Subscript[iim, index]),1},(Subscript[ire, index_]):>{(Subscript[ire, index]),1},(Subscript[ure, index_]):>{(Subscript[ure, index]),Unom},(Subscript[uim, index_]):>{(Subscript[uim, index]),Unom}};
solFindRoot=FindRoot[eqnsTogether,secondParamteterOfFindRoot]//Chop

mf=MatrixForm;
vectfces=First/@eqnsTogether;
vars=unknownVariables;
grad[expr_]:=D[expr,#]&/@vars;
jacobian=grad/@vectfces;
mf@%;
xStart=unknownVariables/.{(Subscript[iim, index_]):>1,(Subscript[ire, index_]):>1,(Subscript[ure, index_]):>Unom,(Subscript[uim, index_]):>Unom};
xNext[xPrevious_]:=xPrevious-Inverse[jacobian/.Thread[vars->xPrevious]].(vectfces/.Thread[vars->xPrevious]);
aux=NestList[xNext,xStart,50][[-1]];
Thread[vars->aux]
solFindRoot

mf@(Chop@jacobian)