@Axel Vogt and @vv 4196

Actually, the main point of my post below of yesterday Aug 10 was to provide the correct answer after eight years of only erroneous ones being stated.

inv_erf(1-10^(-20))=6.60158062235514256151639163241870746485809634382703325836597915047821888510144..

It is usually still of interest of course,  to find other methods, once a problem has been solved. Also I wanted to share an elegant cubic convergent method, variant of Schroeder's series.

      dl:= -f(x)/f'(x);
     dlt := dl-dl^2*f''(x)/(2*f'(x));
     dlt := dl-dlt^2*f''(x)/(2*f'(x));
     x:=x+dlT;

Afterwards, I reran the code with Digits:=249 and it worked fine, converging after six (6) steps of my enhanced Schroeder series method.  The first 80 decimal places are as follows:

         inv_erf(1-10^(-20))=
6.60158062235514256151639163241870746485809634382703325836597915047821888510144...

The iteration step in the 2nd order Taylor expansion is by simple substitution:

     d := -f(x)/f'(x);
     dl := d-d^2*f''(x)/(2*f'(x));
     dl := d-dl^2*f''(x)/(2*f'(x));
     x:=x+d;

If you omit the third line, then you get the Schroeder series.  It is tempting to use the quadratic formula for solving the quadratic equation exactly, but you may run into imaginary roots, as I did in this case. Adding the extra substitution in the third linegreatly reduces the residual in the Taylor expansion.