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: 
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.
PSTricks solution, using package `pst-ode`, command `\pstODEsolve`.

The integral is thought of as an initial value problem and solved (integrated) numerically between 0 and *x* using RKF45 method with automatic stepsize control.

[![enter image description here](https://i.stack.imgur.com/Thqwu.png)](https://i.stack.imgur.com/Thqwu.png)

```
\documentclass[pstricks,border=5pt,12pt]{standalone}
\usepackage{pst-ode,pst-plot}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% erf
% #1: x
% #2: PS variable for result list
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\erf(#1)#2{%
  \pstODEsolve[algebraicAll]{retVal}{ %`retval' takes initial value and erf(#1)
    2/sqrt(Pi)*y[0] %post-multiply result of integration with 2/sqrt(Pi)
  }{0}{#1}{ % integration interval
    2 % two output points (at t=0 and t=#1)
  }{
    0.0 % initial value at t=0
  }{Euler^(-t^2)} %integrand
  %
  %  initialise, if necessary, empty solution list given as arg #2
  \pstVerb{/#2 where {pop}{/#2 {} def} ifelse}
  %  From `retVal', we throw away result at t=0, and append x (arg #1) and erf(x) to our
  %  solution list (arg #2):
  \pstVerb{/#2 [#2 #1 retVal exch pop] cvx def}
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

%solve function integral, 121 plotpoints: -3, -2.95, ..., 2.95, 3
\multido{\nx=-3.00+0.05}{121}{
  \erf(\nx){erf} %result appended to solution list `erf'
}

%plot result
\begin{pspicture}(-1.5,-1)(0.2,6.2)
  \psset{xAxisLabel=$x$,yAxisLabel=$\mathrm{erf}(x)$,xAxisLabelPos={c,-0.3},yAxisLabelPos={-1,c}}
  \begin{psgraph}[axesstyle=frame,Oy=-1,Ox=-3,Dx=1,Dy=0.5](-3,-1)(-3,-1)(3,1){8cm}{6cm}
  \listplot[linecolor=red]{erf}
  \end{psgraph}
\end{pspicture}
\end{document}

```