Let's work it out from first principles, beginning with the hard work of computing a convolution.
As an auxiliary calculation, consider the distribution of $W=X+Y$ where $Y$ has an Exponential distribution with pdf $$f_Y(y) = e^{-y}\,\mathcal{I}(y\gt 0)$$ and $X$ has a Normal$(\mu,\sigma^2)$ distribution with pdf $f_X(x;\mu,\sigma) = \phi((x-\mu)/\sigma)/\sigma$ where $$\phi(z) = \frac{1}{\sqrt{2\pi}}\,e^{-z^2/2}$$ is the standard Normal pdf. The PDF of the sum is the convolution
$$f_W(w;\mu,\sigma) = \int_{-\infty}^\infty f_Y(y) f_X(w-y;\mu,\sigma)\,\mathrm{d}y = \int_0^\infty e^{-y} f_X(w-y;\mu,\sigma)\,\mathrm{d}y.$$
Substituting $\sigma z = w - y - \mu$ expresses this integral as
$$\eqalign{f_W(w;\mu,\sigma) &= e^{\mu-w}\,e^{\sigma^2/2}\int_{-\infty}^{(w-\mu)/\sigma} \phi(z-\sigma)\,\mathrm{d}z \\
&= e^{\mu-w+\sigma^2/2}\, \Phi\left(\frac{w-\mu}{\sigma}-\sigma\right)}\tag{1}$$
where $\Phi$ is the standard normal CDF,
$$\Phi(z) = \int_{-\infty}^z \phi(z)\,\mathrm{d}z.$$
The rest builds on this work and is relatively easy.
An asymmetric Laplace random variable $U$ is based on a mixture of a scaled exponential distribution and the negative of a scaled exponential distribution (potentially with a different scale, thereby making the mixture asymmetric). This mixture is then shifted by a specified amount. The amount of mixing is established to give the Laplace pdf a unique value at its peak--but this is unimportant.
One component of $U$ therefore can be expressed as $$U_+ = \alpha Y + \lambda$$ with a positive scale $\alpha$ and the other component as $$U_- = -\beta Y + \lambda$$ with a positive scale $\beta.$ (I apologize: I worked this out before realizing that my $\alpha$ is $1/\alpha$ in the paper and my $\beta$ is $1/\beta$ in the paper: in the end, after setting $\alpha=\beta,$ this won't matter.)
When we add $X = \sigma Z + \mu$ we obtain two components, of which the first is $$W_+ = U_+ + X = \alpha Y + \lambda + \sigma Z + \mu = \alpha\left(Y + \left[\frac{\sigma}{\alpha} Z + \frac{\lambda + \mu}{\alpha}\right]\right)$$ and the second is similarly written. To obtain its pdf, all we need to do is scale formula $(1)$ by $\alpha,$ giving
$$f_{W_+}(w;\mu,\sigma,\lambda,\alpha) = \frac{1}{\alpha}\,f_W\left(\frac{w}{\alpha};\frac{\lambda+\mu}{\alpha}, \frac{\sigma}{\alpha}\right).\tag{2}$$
Likewise, because
$$W_- = U_- + X = -\beta Y + \lambda + \sigma Z + \mu = -\beta \left(Y + \left[-\frac{\sigma}{\beta } Z + \frac{\lambda + \mu}{\beta }\right]\right)$$
and $-Z$ has the same distribution as $Z$, formula $(1)$ yields
$$\eqalign{f_{W_-}(w;\mu,\sigma,\lambda,\beta) &= \frac{1}{\beta }\,f_W\left(-\frac{w}{\beta };-\frac{\lambda+\mu}{\beta }, \frac{\sigma}{\beta }\right) \\
&= f_{W_+}(-w;-\lambda,\beta,-\mu,\sigma).}\tag{3}$$
The mixture pdf is
$$f_W(w;\mu,\sigma,\lambda,\alpha,\beta,p) = pf_{W_+}(w;\mu,\sigma,\lambda,\alpha) + (1-p) f_{W_-}(w;\mu,\sigma,\lambda,\beta).\tag{4}$$
Comments
For the Laplace-Normal distribution, use $p = \alpha / (\alpha + \beta).$ In your case $\alpha=\beta,$ which evidently "simplifies" $(4)$ a tiny bit--but a quick look at its component formulas $(2)$ and $(3)$ suggests there's not much one can do algebraically to reduce the amount of computation, so why bother?
Each of the components of the final formula $(4),$ as embodied in formulas $(1),$ $(2),$ and $(3)$ can be separately and flexibly implemented and separately tested. This makes for an easier and more reliable software implementation than attempting to combine them all into one monster combination of $\phi$ and $\Phi,$ as done in the referenced paper. As a bonus, important numerical improvements in the calculation can be implemented exactly where they are needed, making the code relatively easy to maintain. As an example, see how $f_W$ is implemented using logarithms (as f.1) in the code below.
Illustration
This plot compares a histogram of one million iid draws from an asymmetric Laplace-Normal distribution with pdf $f_W(w;4,0.5,-3,2,1,2/3)$ to a calculation based directly on formulas $(1) - (4):$
The agreement is a pretty good test.
Code
Here's the R code that generated this simulation and this plot.
n <- 1e6 # Size of simulation
mu <- 4
sigma <- 1/2
alpha <- 2
lambda <- -3
beta <- 1
#
# Generate data.
# set.seed(17)
X <- rnorm(n, mu, sigma)
Y <- ifelse(runif(n, 0, alpha + beta) < alpha, alpha, -beta) * rexp(n) + lambda
W <- X + Y
#
# Plot their histogram.
#
hist(W, freq=FALSE, breaks=200, cex.main=1)
#
# Overplot the PDF.
#
f.1 <- function(w, mu=0, sigma=1) {
exp(mu - w + sigma^2/2 + pnorm((w - mu)/sigma - sigma, log=TRUE))
}
f.plus <- function(w, mu=0, sigma=1, lambda=0, alpha=1) {
f.1(w / alpha, (lambda + mu) / alpha, sigma / alpha) / alpha
}
f.minus <- function(w, mu=0, sigma=1, lambda=0, beta=1) {
f.plus(-w, -mu, sigma, -lambda, beta)
}
f <- function(w, mu=0, sigma=1, lambda=0, alpha=1, beta=1, p=1/2) {
p * f.plus(w, mu, sigma, lambda, alpha) + (1-p) * f.minus(w, mu, sigma, lambda, beta)
}
f.asymmetric <- function(y, mu=0, sigma=1, lambda=0, alpha=1, beta=1) {
f(y, mu, sigma, lambda, alpha, beta, alpha / (alpha + beta))
}
curve(f.asymmetric(x, mu, sigma, lambda, alpha, beta), add=TRUE, lwd=2, col="Red")
