This integral can be found in *D. B. Owen (1980) A table of normal integrals, Communications in Statistics - Simulation and Computation, 9:4, 389-419*: [![enter image description here][1]][1] `BvN` denotes the bivariate normal probability function. Check in R: > a <- 2 > b <- 3 > w <- 5 > f <- function(x) dnorm(x)*pnorm(a+b*x) > integrate(f, lower=-Inf, upper=w) 0.7364551 with absolute error < 1.3e-06 > > rho <- -b/sqrt(1+b^2) > Sigma <- cbind(c(1,rho),c(rho,1)) > mvtnorm::pmvnorm(upper=c(a/sqrt(1+b^2), w), sigma=Sigma) [1] 0.7364551 attr(,"error") [1] 1e-15 attr(,"msg") [1] "Normal Completion" Alternatively, you can express this integral with the Owen $T$-function: > library(OwenQ) > 1/2*(pnorm(a/sqrt(1+b^2)) + pnorm(w) - 2*OwenT(w, (b*w+a)/w) - 2*OwenT(-a/sqrt(1+b^2), (a*b+w*(1+b^2))/a) - (a <= 0)) [1] 0.7364551 [1]: https://i.sstatic.net/fIOSQ.png