Instead of a PDE solution we can examine the dual problem and integrate by various versions of Romberg integraton as a 2d loop. This works showing that the integrand is well-behaved.
Looking at 26.3.3 we can see that the inner integral can be directly written in terms of C++11 erf(x). Then in the other direction we can use Romberg/Midpoint to get accuracy. In this case we don't need to build a matrix. And C++11 takes care of all edge conditions in erf(x) and exp(x). Just in case you can do a Kahan summation.
http://people.math.sfu.ca/~cbm/aands/page_936.htm
Some results starting with a modest N = 10 gives this (the result are 7 digits accuracy without too much effort). N = 50,100, 400 are stress tests.
So we can get accurate results directly in C++11 without external libraries. It's another option.
a, b, rho: 1.93185,1.78999, -0.411002
*Genz West : 0.9366216225051415
*Tanh 2d Extrap/Adaptive : 0.9366219402518232
*Midpoint 2d : 0.9366238179854843
*Tanh 2d 100X100 : 0.9367276920000295
*A&S 26.3.3 : 0.9366704332207069
*Genz QuantLib 1.8 : 0.9366216225051415
*A&S 26.3.3 Extrap (N=10) : 0.9366210429626975
*A&S 26.3.3 Extrap (N=50) : 0.9366216224739409
*A&S 26.3.3 Extrap (N=100) : 0.9366216225046546
*A&S 26.3.3 Extrap (N=400) : 0.9366216225051428
a, b, rho: 1.337236842557504,1.264162992646285, 0.486393246573539
*Genz West : 0.8363779098801795
*Tanh 2d Extrap/Adaptive : 0.8363726930947498
*Midpoint 2d : 0.836377419248175
*Tanh 2d 100X100 : 0.8365230433697367
*A&S 26.3.3 : 0.8364483519381961
*Genz QuantLib 1.8 : 0.8363779098801796
*A&S 26.3.3 Extrap (N=10) : 0.8363776712097363
*A&S 26.3.3 Extrap (N=50) : 0.836377909871765
*A&S 26.3.3 Extrap (N=100) : 0.8363779098800628
*A&S 26.3.3 Extrap (N=400) : 0.8363779098801919
a, b, rho: 1.955091230607987,0.941464704072475, -0.4825543785261079
*Genz West : 0.80173352113734
*Tanh 2d Extrap/Adaptive : 0.8017309528547933
*Midpoint 2d : 0.8017282455937149
*Tanh 2d 100X100 : 0.8018595575746997
*A&S 26.3.3 : 0.8017798723119439
*Genz QuantLib 1.8 : 0.80173352113734
*A&S 26.3.3 Extrap (N=10) : 0.8017329230670884
*A&S 26.3.3 Extrap (N=50) : 0.8017335211054371
*A&S 26.3.3 Extrap (N=100) : 0.8017335211368406
*A&S 26.3.3 Extrap (N=400) : 0.8017335211373599
a, b, rho: 1.860920640238858,0.3513236324632962, -0.09898276912896131
*Genz West : 0.614801617049458
*Tanh 2d Extrap/Adaptive : 0.6147996574002191
*Midpoint 2d : 0.6147982757717783
*Tanh 2d 100X100 : 0.6148746733453729
*A&S 26.3.3 : 0.6148381569839125
*Genz QuantLib 1.8 : 0.6148016170494581
*A&S 26.3.3 Extrap (N=10) : 0.6148012461785963
*A&S 26.3.3 Extrap (N=50) : 0.6148016170304943
*A&S 26.3.3 Extrap (N=100) : 0.6148016170491677