In a comment to a Mapleprimes thread, Jacques mentioned an old suggestion of Kahan's that numerical computations should return an estimate of conditioning alongside a result.

I mentioned in this comment an approach for numerical estimation of (all) roots of a univariate polynomial with real or complex numeric coeffficients that is based upon computing eigenvalues of a companion matrix. Here below is some rough code to inplement that idea, but which also returns condition number estimates associated with the eigenvalues.

I include an example of the badly conditioned Wilkinson's polynomial. It is possible that better results could be obtained by using a Lagrange basis representation of that polynomial, but I didn't try to figure out how that would work in an analogous way. The standard Maple utility, fsolve, has no problem with this example.

froots := proc(
  p::polynom
, {a::complex(extended_numeric):=-infinity-infinity*I}
, {b::complex(extended_numeric):=infinity+infinity*I}
, {allowcomplex::truefalse:=false}
, $
)
local A,boundar,boundai,boundbr,boundbi,cond,C,dtype,eps,N,sol,var;
var := indets(p)[1];
if has(p,I) then
  dtype:='complex'('float');
else
  dtype:='float';
end if;
if not allowcomplex then
  if not type(a,extended_numeric) then
    (boundar,boundai):=Re(a),0;
  else
    (boundar,boundai):=Re(a),Im(a);
  end if;
  if not type(b,extended_numeric) then
    (boundbr,boundbi):=Re(b),0;
  else
    (boundbr,boundbi):=Re(b),Im(b);
  end if;
else
  (boundar,boundai):=Re(a),Im(a);
  (boundbr,boundbi):=Re(b),Im(b);
end if;
A := LinearAlgebra:-LA_Main:-CompanionMatrix(
            evalf(p), var, 'compact'=false,
            'outputoptions'=['datatype'=dtype]);
if nops([A])=1 and coeff(p,var,degree(p,var)) = 1 then
  N := op([1,1],A);
  C :=  Matrix(N,N,'storage'='empty','shape'='identity','datatype'=dtype);
else
  (A,C):=A;
  N := op([1,1],A);
end if;
eps:=Float(5,-Digits);
(sol,cond) := LinearAlgebra:-LA_Main:-EigenConditionNumbers(
         A, C,
         'balance'=both,
         ':-output'=[':-values',':-conditionvalues'],
         'outputoptions[conditionvalues]'=[],
         'outputoptions[conditionvectors]'=[],
         'outputoptions[values]'=['datatype'='complex'('float')],
         'outputoptions[vectors]'=[]);
sol := select(x -> evalb(Re(x) <= boundbr and Im(x) <= boundbi
                         and boundar <= Re(x) and boundai <= Im(x)),
              evalf(convert(sol,list)));
seq(`if`(abs(Im(sol[i])) < eps,[Re(sol[i]),evalf(cond[i])],
         [sol[i],evalf(cond[i])]), i=1..nops(sol));
end proc:
 
# Wilkinson's polynomial.
y := expand(mul((x-i),i=1..20)):
                                                                                
st:=time():fsolve(y,x);time()-st;
st:=time():froots(y);time()-st;
                                                                                
st:=time():fsolve(y,x,complex);time()-st;
st:=time():froots(y,allowcomplex=true);time()-st;

I could alter the section of code which produces the companion matrix to be more efficient. But I have a hope to figure out better how to use CompanionMatrix here with a Lagrange basis. I should probably look harder at the code in RootFinding[BivariatePolynomial].

acer