Here it is as an `R` function `f`: `
f <- function(x, lambda) {
pnorm(sqrt(x+1/2), sqrt(lambda), 1/2) - pnorm(sqrt(pmax(0, x-1/2)), sqrt(lambda), 1/2)
};
lambda <- 3;
x <- 0:(3*lambda);
y <- dpois(x, lambda);
i <- y > 1e-6;
plot(x[i], y[i], pch=21, bg="#e02020");
lines(x[i], f(x[i], lambda), type="h", add=TRUE)`