Question: Is there a way to easily compute [the erf function](https://en.wikipedia.org/wiki/Normal_distribution#Cumulative_distribution) (or the cumulative distribution function of the normal law) in LaTeX?

Currently, I use `pgf` to make computation, but I did not find a way to compute erf using `pgf`.

I would be happy to use any package that is available to compute erf, or any custom solution to compute that function.

Answer: 
For precise values, I recommend externalizing the calculation, here `gnuplot` is used.

Code (needs `--shell-escape` enabled)
-------------------------------------

```
\documentclass{article}
\usepackage{amsmath,pgfmath,pgffor}
\makeatletter
\def\qrr@split@result#1 #2\@qrr@split@result{\edef\erfInput{#1}\edef\erfResult{#2}}
\newcommand*{\gnuplotErf}[2][\jobname.eval]{%
    \immediate\write18{gnuplot -e "set print '#1'; print #2, erf(#2);"}%
    \everyeof{\noexpand}
    \edef\qrr@temp{\@@input #1 }%
    \expandafter\qrr@split@result\qrr@temp\@qrr@split@result
}
\makeatother
\DeclareMathOperator{\erf}{erf}
\begin{document}
\foreach \x in {-50,...,50}{%
\pgfmathparse{\x/50}%
\gnuplotErf{\x/50.}%
$ x = \pgfmathresult = \erfInput, \erf(x) = \erfResult$\par
}
\end{document}

```

Output
------

![enter image description here](https://i.stack.imgur.com/hThhI.png)
Using the [same approximation idea as cjorssen](https://tex.stackexchange.com/a/107439/13450) (I tried the Taylor series as Qrrbrbirlbel suggested but it's pretty hopeless to get a decent approximation this way) I rewrote the function without using low-level PGF. Because we have so many 2D plots here already, I'll just use my 3D plot that I had already.

```
\documentclass{standalone}
\usepackage{pgfplots}
\usepackage{tikz}
\pgfplotsset{
colormap={bluewhite}{ color(0cm)=(rgb:red,18;green,64;blue,171); color(1cm)=(white)}
}
\begin{document}
\begin{tikzpicture}[
    declare function={erf(\x)=%
      (1+(e^(-(\x*\x))*(-265.057+abs(\x)*(-135.065+abs(\x)%
      *(-59.646+(-6.84727-0.777889*abs(\x))*abs(\x)))))%
      /(3.05259+abs(\x))^5)*(\x>0?1:-1);},
    declare function={erf2(\x,\y)=erf(\x)+erf(\y);}
]
\begin{axis}[
    small,
    colormap name=bluewhite,
    width=\textwidth,
    enlargelimits=false,
    grid=major,
    domain=-3:3,
    y domain=-3:3,
    samples=33,
    unit vector ratio*=1 1 1,
    view={20}{20},
    colorbar,
    colorbar style={
        at={(1,-.15)},
        anchor=south west,
        height=0.25*\pgfkeysvalueof{/pgfplots/parent axis height},
    }
]
\addplot3 [surf,shader=faceted] {erf2(x,y)};
\end{axis}
\end{tikzpicture}
\end{document}

```

![3D plot of erf(x)+erf(y)](https://i.stack.imgur.com/SfCxc.png)

The approximation has a maximum error of 1.5·10-7 ([source](http://en.wikipedia.org/wiki/Error_function#Approximation_with_elementary_functions)).

Thanks to [Jake](https://tex.stackexchange.com/users/2552/jake) for [spotting and fixing](http://chat.stackexchange.com/transcript/message/9964403#9964403) the wrong syntax I first had in this code.