By the way: using a 1F1 representation for erf, the values agree with Wolfram | Alpha. But the hypergeometric code has to switch to bigfloats to do this:

(%i11) my_erf(x) := nfloat(2*x*hypergeometric([1/2],[3/2],-x^2)/sqrt(%pi),[]);

(%o11) my_erf(x):=nfloat((2*x*hypergeometric([1/2],[3/2],-x^2))/sqrt(%pi),[])

(%i12) load(hypergeometric)$

(%i14) my_erf(1.0+5*%i);

(%o14) -2.783777029220896b9*%i-1.078693116198541b9