Just a few remarks about the code in "InvErf_1.vi":
1) Because of the Fundamental Theorem of Calculus, we can compute an explicit derivative here, and there is no need for resorting to finite difference ratios.
In other words, because of the FTC, the derivative of erf(x) is known to be
[2 / sqrt(pi)] exp(-x2)
2) It is very easy to create [admittedly pathological] functions which cause an infinite loop in Newton-Raphson-Simpson. For instance, you can create a scenario with "Twin Peaks", and a zero between them, but with slopes near each peak which point off towards the other peak, so that if you make the mistake of starting immediately beneath one of those special points, then Newton-Raphson-Simpson will just bounce you back and forth forever between the two peaks, and you'll never work your way down to the zero in the valley.
Bottom line being that because Newton-Raphson-Simpson can give you infinite loops [even in the presence of a known zero], it's imperative that you put some sort of upper bound on the number of iterations in your "while" loop so that your program won't spin out of control if it comes into contact with this sort of pathological behavior.
3) This point is very subtle, but bear with me: Because of the [ABSOLUTELY MADDENING!!!] nature of floating-point numbers [as opposed to fixed-point numbers], THE GRANULARITY OF FLOATING-POINT NUMBERS IS NON-CONSTANT!!!
As you move up and down the "Real Axis", the granularity of floating-point numbers changes rather drastically - down around zero, a floating-point implementation will have fairly good [small, tight] granularity, but as you move out towards positive infinity, the granularity quickly bulges way out to the point that it becomes useless for any serious numerical computations.
In your algorithm, you've got both a constant "accuracy" and a constant "derivative increment", which shouldn't be much of a problem in the general vicinity of regions of known granularity, but if you find yourself pushed out into regions of possibly lousy [large, loose] granularity, then it's easy to create examples where an algorithm like yours could [at least in theory] find itself in an infinite loop because the granularity of the fixed-point numbers is locally much larger than the constant increments against which you are comparing the differences of such fixed-point numbers. For instance, it's easy to fiddle with the "Twin Peaks" example above so that an approach to Newton-Raphson-Simpson like yours will converge in regions of tight floating-point granularity, but when translated out to regions of loose floating-point granularity, will find itself locked in an infinite loop.
In the case of Newton-Raphson-Simpson as it's used in the Inverse Function Theorem, this is particularly infuriating because you may not know a priori where the inverse came from - it could be that the function value lives in a region with a very tight granularity of floating-point numbers, but that the inverse lives way off in some faraway region with very loose granularity of floating-point numbers.
Anyway, if you ever get really serious about numerical analysis, then you need to think about compensating for variable granularity in an algorithm like this [or proving a priori that you are necessarily confined to regions of tight granularity, so that any potential problems with granularity can be ignored].
4) Finally, while erf(x) doesn't have the kinds of "Twin Peaks" pathologies that I used above [erf(x) being monotonically increasing], it is nevertheless the case that in both the positive and negative directions, erf(x) VERY QUICKLY moves up towards its horizontal asymptotes.
Horizontal asymptotes, in turn, mean that the slopes f'(xn) are heading down towards zero [as above, by the FTC, the slopes are moving towards zero at the order of exp(-x2)], and at some point, these slopes could very well be rounded down to zero [and with non-constant granularity of the floating-point numbers, it's a real pain-in-the-neck to know exactly when the metal will decide to round the slopes down to zero]
And once you get a zero derivative in Newton-Raphson-Simpson & the Inverse Function Theorem, then you've got one heckuva mess on your hands.
So if you find yourself trying to compute the inverse error function for error function values down towards -1, or up towards +1, then you need to be very careful that the metal doesn't round your derivatives down to zero and completely ruin your algorithm.
