I coded simpler (and *MUCH* wronger) versions of this general code countless times…

Thank you Steven for the marvellous work! register it ASAP, someone much wiser than me said: “its usefulness for practical computations can hardly be overestimated.”

EDIT 1:

It can be useful to protect tricky computations in AD:

```
using Richardson
using ForwardDiff: Dual
f(x) = cos(x + sqrt(x) * x)
A,_=extrapolate(Dual(0.1,0.0),x0=Dual(0.0,1.0)) do x
@show x
f(x)
end
#=
x = Dual{Nothing}(0.1,2.0)
x = Dual{Nothing}(0.0125,1.125)
x = Dual{Nothing}(0.0015625,1.015625)
x = Dual{Nothing}(0.0001953125,1.001953125)
x = Dual{Nothing}(2.44140625e-5,1.000244140625)
x = Dual{Nothing}(3.0517578125e-6,1.000030517578125)
=#
B=f( Dual(0.0,1.0) )
@show A # Dual{Nothing}(0.9999999999995401,1.740011327412078e-8)
@show B # Dual{Nothing}(1.0,NaN))
```

The following is now fixed/supported!

```
#EDIT 2: via super ugly, surely broken approach it seems to work for vector valued functions, is there a way to truly support vector functions?
#```julia
using LinearAlgebra: norm
import Base.abs
import Base.isless
Base.abs(x::AbstractVector) = norm(x)
Base.isless(y::Number,x::AbstractVector) = y<norm(x)
extrapolate(-0.1,x0=0) do x
@show x
[x/sqrt(-x),sin(x)/x]
end
# ([-4.76091483644921e-9, 1.0], 8.704985825570316e-9)
#```
#EDIT 3: limits in R^N can be easily supported i think
# ```julia
using Richardson
@inline function Richardson.extrapolate(f,h::AbstractArray{T};x0::AbstractArray{T},args...) where T
@boundscheck size(h)==size(x0)
@inline g(t)=@inbounds x0.+t.*h
extrapolate(T(1),x0=T(0);args...) do t
f(g(t))
end
end
extrapolate([0.1,-0.2],x0=[0.0,1.0]) do x
@show x
sin(x[1])/x[1]*log(x[2])/(x[2]-1)
end
#=
x = [0.1, 0.8]
x = [0.0125, 0.975]
x = [0.0015625, 0.996875]
x = [0.0001953125, 0.999609375]
x = [2.44140625e-5, 0.999951171875]
=#
#(1.000000000000058, 1.4038002982275088e-9)
#```
#EDIT 4: for some functions it can go into endless loops (valid up to commit: https://github.com/JuliaMath/Richardson.jl/commit/c828d05a0b07fde094092496eecdea0e1c4588c9)
#```julia
using Richardson
extrapolate(0.2,x0=1.0) do x
@show x # this is there to slow it down just enough to be able to abort the computation easily
log(x)/sqrt(x-1)
end
#```
```