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*)

SeedRandom[1]
n=3;
{a,b,c}=randomMatrix[n]&/@Range[n];
{p1,p2,p3}=-3Random[]&/@Range[n];
f1[xvect_]:=Plus@@(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])];
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]
Solve[origEqns]
(*for linear equations one step is enough*)

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]