Numerical Instability with this cumulate I am trying to compute this cumulate. However, with very high or very low values of x, the computation is affected by numerical errors.
$$p=1-\frac{2 u}{\sqrt{2 \pi } l \sqrt{t} e^{\frac{(l t+u x)^2}{2 t u^2}} \left(\text{erf}\left(\frac{l t+u x}{\sqrt{2} \sqrt{t} u}\right)-\text{erf}\left(\frac{l t+u (x-t u)}{\sqrt{2} \sqrt{t} u}\right)\right)+2 u}$$
which, for t=1; u=0.5 and l=2 gives this shape which is clearly affected by numerical errors on both sides:

I thinks the problem is due to the really high number in the exponential, which is then multiplied by the really small number of the to erf functions. I tried different solutions (mostly involving expanding the exponential) but none worked.
The computation are done in MATLAB
Any help is appreciated.
 A: Try writing the denominator in terms of the complementary error function $\mathrm{erfc}(x) = 1 - \mathrm{erf}(x)$. Doing so should give you more accuracy for the larger positive values of $x$. The negative values of $x$, however, are still a problem.
Here is a better solution. The function of $x$ in the denominator can be written as
$$
f(x) = \frac{2}{\sqrt{\pi}}\int_a^b \exp(-(\tau + b)(\tau - b))\,d\tau,
$$
where $a = \alpha + \beta(x - tu)$, $b = \alpha + \beta x$, $\alpha = \ell t /(u\sqrt{2t})$, and $\beta = 1/\sqrt{2t}$. It follows that
$$
p = 1 - \frac{2u}{\ell\sqrt{2\pi t}f(x) + 2u}
$$
For each value of $x$, use a quadrature routine with a small convergence tolerance to evaluate $f$. Doing so, you can accurately compute $p(x)$ over a wider range of values. Clearly this approach is more computationally demanding then using $\mathrm{erf}$ or $\mathrm{erfc}$, however, it should be more robust to the input value $x$. In practice, you should be able to calculate $p$ for ''not too negative'' values of $x$ using the solution mentioned above and then resort to this approach for negative values of $x$ past some cutoff. Here is a plot I made in Matlab using their $\texttt{quad}$ function to evaluate $f$ with a tolerance of $10^{-12}$, $x\in[-15,20]$, and your stated parameter values.

