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: 
Based on [this answer](https://stackoverflow.com/a/457805/1195799).

```
\documentclass{standalone}
\usepackage{tikz}
\makeatletter
\pgfmathdeclarefunction{erf}{1}{%
  \begingroup
    \pgfmathparse{#1 > 0 ? 1 : -1}%
    \edef\sign{\pgfmathresult}%
    \pgfmathparse{abs(#1)}%
    \edef\x{\pgfmathresult}%
    \pgfmathparse{1/(1+0.3275911*\x)}%
    \edef\t{\pgfmathresult}%
    \pgfmathparse{%
      1 - (((((1.061405429*\t -1.453152027)*\t) + 1.421413741)*\t 
      -0.284496736)*\t + 0.254829592)*\t*exp(-(\x*\x))}%
    \edef\y{\pgfmathresult}%
    \pgfmathparse{(\sign)*\y}%
    \pgfmath@smuggleone\pgfmathresult%
  \endgroup
}
\makeatother
\begin{document}
\begin{tikzpicture}[yscale = 3]
  \draw[very thick,->] (-5,0) -- node[at end,below] {$x$}(5,0);
  \draw[very thick,->] (0,-1) -- node[below left] {$0$} node[at end,
  left] {$erf(x)$} (0,1);
  \draw[red,thick] plot[domain=-5:5,samples=200] (\x,{erf(\x)});
\end{tikzpicture}
\end{document}

```

![enter image description here](https://i.stack.imgur.com/VkEZI.png)
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}

```