There have been a few posts on mapleprimes about numerically solving systems of procedures. The latest one, up until now, was this.

Here's some code to implement the method. Since the algorithm is basically very simple, I've added a few bells and whistles as optional arguments.

The essence of it is as follows. The number of procedures must match the number of parameters of each and every procedure. It does maxtries attempts at choosing a random point, and then does at most maxiter iterations. A solution is only accepted if the norm of the last change in vector (point) x is less than xtol, and if the forward error norm(F(x)) is less than ftol. The jacobian of F may be supplied optionally as a Matrix of procedures, or a method for computing the jacobian may be supplied. The methods are fdiff which only uses Maple's numerical differentiation routine fdiff, or hybrid which attempts symbolic differentiation via Maple's D[] operator and then falls back to fdiff via the nifty evalf@D equivalence.

I have not profiled it, though I have attempted to set up the jacobian construction so that it re-uses the same rtable for each instantiation. I have not gone really low-level with the external-calling to numeric solvers. I have not tried to leverage fast hardware double-precision solutions as jumping-off points for extended precision polishing (eg. see here).

I've added some examples, after the code below. Feel free to comment or mention bugs. I'd really like a better name for it.

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}):=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:

Here are a few examples, which I have not set up as verified checks or tests.

kernelopts(printbytes=false):
kernelopts(gcfreq=10^8);
 
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;
 
NRprocs([f,g,(a,b)->a+b]);
NRprocs([(a,b,c)->a*b*c]);
 
fl:=[x->cos(x)-x]:
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]:
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])):
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:=[(x,y,z)->x*z+2.1*y^2-sqrt(z),(x,y,z)->z+1.1*x*y,(x,y,z)->x^3+1.1*y*z]:
gc():
st,ba,bu:=time(),kernelopts(bytesalloc),kernelopts(bytesused):
for i from 1 to 5 do
  X0:=LinearAlgebra:-RandomVector(3,generator=-1.0..1.0,outputoptions=[datatype=
float]);
  NRprocs(fl,initialpoint=X0);
end do:
time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;
 
fl:=[(x,y,z)->x*z+2.2*y^2-sqrt(z),(x,y,z)->z+1.2*x*y,(x,y,z)->x^3+1.2*y*z]:
gc():
st,ba,bu:=time(),kernelopts(bytesalloc),kernelopts(bytesused):
jfl:=Matrix(3,3,(i,j)->D[j](fl[i])):
for i from 1 to 5 do
  X0:=LinearAlgebra:-RandomVector(3,generator=-1.0..1.0,outputoptions=[datatype=
float]);
  NRprocs(fl,jacobian=jfl,initialpoint=X0);
end do:
time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;
 
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]:
gc():
st,ba,bu:=time(),kernelopts(bytesalloc),kernelopts(bytesused):
jfl:=Matrix(3,3,(i,j)->D[j](fl[i])):
time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;
sol:=NRprocs(fl,jacobian=jfl,maxtries=50);
time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;
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]:
gc():
st,ba,bu:=time(),kernelopts(bytesalloc),kernelopts(bytesused):
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):
time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;
sol:=NRprocs(cfl,jacobian=cjfl);
time()-st,kernelopts(bytesalloc)-ba,kernelopts(bytesused)-bu;
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);
NRprocs([f,g],method=fdiff);

acer