restart:

NRprocs:=proc(

  funlist::{list(procedure),Vector(procedure)}

, {allowcomplex::truefalse:=false}

, {initialpoint::Vector:=NoUserValue}

, {maxtries::posint:=NoUserValue}

, {maxiter::posint:=30}

, {xtol::realcons:=5.0*10^(-Digits)}

, {ftol::realcons:=1.0*10^(-trunc(2*Digits/3))}

, {method::{identical(hybrid,fdiff)}:=':-hybrid'}

, {jacobian::{Matrix({operator,procedure}),identical(NoUserValue)}:=':-NoUserValue'}

, $

)

local a, b, currXseq, dummy, dummyseq, F, i, ii, J, j, jj,

      MaxTries, numvars, oldX, thisJ, thistry, varnumset, X;

  varnumset := {seq(`if`(not(member('call_external',{op(3,eval(F))})),
                         nops([op(1,eval(F))]),NULL), F in funlist)};

  if nops(varnumset) > 1 then

    error "expecting procedure all taking same number of arguments";

  elif nops(varnumset) = 1 then

    numvars := varnumset[1];

  else

    numvars := nops(funlist);

  end if;

  if numvars <> nops(funlist) then

    error "number of procedures %1 does not match number of their parameters %2"

, nops(funlist), numvars;

  end if;

  if initialpoint <> 'NoUserValue' then

    oldX := Vector(numvars,(i)->evalf(initialpoint[i]),':-storage'=':-rectangular',

                   ':-datatype'=`if`(allowcomplex,'complex'(':-float'),':-float'));

    if maxtries = ':-NoUserValue' then

      MaxTries := 1;

    elif maxtries > 1 then

      error "maxtries must be 1 when initialpoint is supplied";

    else

      MaxTries := maxtries;

    end if;

  else

    if maxtries = ':-NoUserValue' then

      MaxTries := 20;

    else

      MaxTries := maxtries;

    end if;

  end if;

  if jacobian = ':-NoUserValue' then

    # This is only built once, so there should be a way to make

    # it take a little more time and make a J that evaluates quicker.

    # Using evalhf is just one possibility.

    dummyseq:=seq(dummy[i],i=1..numvars);

    if method=':-hybrid' then

      J:=Matrix(nops(funlist),numvars,

              (i,j)->unapply('evalf'(D[j](funlist[i])(dummyseq)),[dummyseq]));

    else # method=fdiff

      J:=Matrix(nops(funlist),numvars,

              (i,j)->unapply(fdiff(funlist[i],[j],[dummyseq]),[dummyseq]));

    end if;

  else

    J:=jacobian;

  end if;

  thisJ := Matrix(nops(funlist),numvars,

                  ':-datatype'=`if`(allowcomplex,'complex'(':-float'),':-float'));

  X := Vector(numvars,':-datatype'=`if`(allowcomplex,'complex'(':-float'),':-float'));

  userinfo(1,'NRprocs',`maximal number of tries is `, MaxTries);

  for thistry from 1 to MaxTries do

    if thistry > 1 or initialpoint = ':-NoUserValue' then

      if allowcomplex = true then

        oldX := LinearAlgebra:-RandomVector(numvars,

                              ':-density'=1,':-generator'=-1.0-1.0*I..1.0+1.0*I,

                              ':-outputoptions'=[':-datatype'='complex'(':-float')]);

      else

        oldX := LinearAlgebra:-RandomVector(numvars,

                              ':-density'=1,':-generator'=-1.0..1.0,

                              ':-outputoptions'=[':-datatype'=':-float']);

      end if;

    end if;

    userinfo(1,'NRprocs',`initial point : `, convert(oldX,list));

    userinfo(1,'NRprocs',`maximal number of iterations is `, maxiter);

    try

      for ii from 1 to maxiter do

        currXseq:=seq(oldX[jj],jj=1..numvars);

        for i from 1 to nops(funlist) do

          for j from 1 to numvars do

            thisJ[i,j] := J[i,j](currXseq);

          end do;

        end do;

        # copy oldX into X, so that increment can be added to X inplace.

        ArrayTools:-Copy(numvars,oldX,0,1,X,0,1);

        LinearAlgebra:-LA_Main:-VectorAdd(

               X, LinearAlgebra:-LA_Main:-LinearSolve(

                   thisJ,

                   Vector(nops(funlist),

                     (i)->evalf(funlist[i](currXseq))),

                   ':-method'=':-none',':-free'=':-NoUserValue',':-conjugate'=true,

                   ':-inplace'=false,':-outputoptions'=[],':-methodoptions'=[]),

               1,-1,':-inplace'=true,':-outputoptions'=[]);

        userinfo(2,'NRprocs',`at iteration`,ii, convert(X,list));

        if LinearAlgebra:-LA_Main:-Norm(

             # oldX is no longer needed, so can overwrite it inplace.

             LinearAlgebra:-LA_Main:-VectorAdd(

               oldX,X,-1,1,':-inplace'=true,':-outputoptions'=[]),

             infinity,':-conjugate'=true) < xtol

         and max(seq(abs(F(seq(X[jj],jj=1..numvars))),F in funlist)) < ftol then

          return X;

        end if;

        ArrayTools:-Copy(numvars,X,0,1,oldX,0,1);

      end do:

    catch "singular matrix","inconsistent system","unable to store":

      ii := ii - 1;

      next;

    end try;

  end do;

  return 'procname'(args);

end proc:

f:=proc(x,y) x^3-sin(y); end proc:

g:=proc(x,y)

local sol;

  sol:=fsolve(q^5+x*q^2+y*q=0,q,complex):

  Re(sol[1]+1);

end proc:

fl:=[f,g]:

NRprocs(fl);

if type(%,Vector) then

  max(seq(abs(fl[i](seq(%[j],j=1..nops(fl)))),i=1..nops(fl)));

end if;

Digits:=32:

NRprocs(fl);

if type(%,Vector) then

  max(seq(abs(fl[i](seq(%[j],j=1..nops(fl)))),i=1..nops(fl)));

end if;

Digits:=10:

infolevel[NRprocs]:=1:

sol := NRprocs([f,g],initialpoint=<5.0, 6.0>);

if type(%,Vector) then

  max(seq(abs(fl[i](seq(%[j],j=1..nops(fl)))),i=1..nops(fl)));

end if;

infolevel[NRprocs]:=0:

infolevel[NRprocs]:=2:

sol := NRprocs(fl,initialpoint=<-1, -2>);

if type(%,Vector) then

  max(seq(abs(fl[i](seq(%[j],j=1..nops(fl)))),i=1..nops(fl)));

end if;

infolevel[NRprocs]:=0:

sol := NRprocs(fl);

if type(%,Vector) then

  max(seq(abs(fl[i](seq(%[j],j=1..nops(fl)))),i=1..nops(fl)));

end if;

sol := NRprocs(fl,ftol=1e-8);

if type(%,Vector) then

  max(seq(abs(fl[i](seq(%[j],j=1..nops(fl)))),i=1..nops(fl)));

end if;

sol := NRprocs(fl,initialpoint=<-I,I>,allowcomplex=true);

if type(%,Vector) then

  max(seq(abs(fl[i](seq(%[j],j=1..nops(fl)))),i=1..nops(fl)));

end if;

sol := NRprocs(fl,allowcomplex=true);

if type(%,Vector) then

  max(seq(abs(fl[i](seq(%[j],j=1..nops(fl)))),i=1..nops(fl)));

end if;

# deliberate errors
NRprocs([f,g,(a,b)->a+b]);

NRprocs([(a,b,c)->a*b*c]);

fl:=[x->cos(x)-x]:

CodeTools:-Usage( NRprocs(fl) );

if type(%,Vector) then

  max(seq(abs(fl[i](seq(%[j],j=1..nops(fl)))),i=1..nops(fl)));

end if;

fl:=[(x,y,z)->x*z+2*y^2-sqrt(z),(x,y,z)->z+x*y,(x,y,z)->x^3+y*z]:

CodeTools:-Usage( NRprocs(fl) );

if type(%,Vector) then

  max(seq(abs(fl[i](seq(%[j],j=1..nops(fl)))),i=1..nops(fl)));

end if;

jfl:=Matrix(3,3,(i,j)->D[j](fl[i])):

CodeTools:-Usage( NRprocs(fl,jacobian=jfl) );

if type(%,Vector) then

  max(seq(abs(fl[i](seq(%[j],j=1..nops(fl)))),i=1..nops(fl)));

end if;

fl:=[proc(x::float,y::float,z::float) x*z+1.9*y^2-z^2; end proc,

     proc(x::float,y::float,z::float) z+0.9*x*y; end proc,

     proc(x::float,y::float,z::float) x^3+0.9*y*z; end proc]:

st,ba,bu:=time(),kernelopts(bytesalloc),kernelopts(bytesused):

jfl:=Matrix(3,3,(i,j)->D[j](fl[i])):

sol:=CodeTools:-Usage( NRprocs(fl,jacobian=jfl,maxtries=50) );

if type(sol,Vector) then

  max(seq(abs(fl[i](seq(sol[j],j=1..nops(fl)))),i=1..nops(fl)));

end if;

fl:=[proc(x::float,y::float,z::float) x*z+1.8*y^2-z^2; end proc,

     proc(x::float,y::float,z::float) z+0.8*x*y; end proc,

     proc(x::float,y::float,z::float) x^3+0.8*y*z; end proc]:

jfl:=Matrix(3,3,(i,j)->D[j](fl[i])):

cfl:=[seq(Compiler:-Compile(fl[i]),i=1..nops(fl))]:

cjfl:=map(t->`if`(type(t,procedure),Compiler:-Compile(t),t),jfl):

sol:=CodeTools:-Usage( NRprocs(cfl,jacobian=cjfl) );

if type(sol,Vector) then

  max(seq(abs(fl[i](seq(sol[j],j=1..nops(fl)))),i=1..nops(fl)));

end if;

NRprocs([f,g],method=hybrid);

[f, g](seq(a, a=%)); # check whether it's a root

NRprocs([f,g],method=fdiff);

[f, g](seq(a, a=%)); # check whether it's a root