A simple way of computing error function is to use Kummer's equation of the form of,
$$ M(a,b,z)=\sum_{s=0}^{\infty} \frac{(a)_s}{(b)_s s!}z^s=1+\frac{a}{b}z+\frac{a(a+1)}{b(b+1)2!}z^2+...,\quad z \in \mathbb{C} $$
and
$$ M(a,b,z)=\sum_{s=0}^{\infty}\frac{(a)_s}{\Gamma{(b+s)}s!}z^s $$
and
$$ \text{erf}(z)=\frac{2z}{\sqrt{\pi}}M(.5,1.5,-z^2)=\frac{2z}{\sqrt{\pi}}e^{-z^2}M(1,1.5,z^2) $$
Here is the R code,
f<-function(a,b,z,maxt=5){ s=1:maxt a=a+s b=b+s ss=1 for(i in s){ mt=prod(a[1:i]/b[1:i]) ss=ss+mt *z^(2*i)/factorial(i) } ss=2*z/sqrt(pi)*exp(-z^2)*ss return(ss) }
here is the plot (for real z) 