I am trying to compute a 2 Dimensional numerical integration with the help of Boost. I've been reading the documentation for a couple of days and my CPP language knowledge is just up to legacy CPP level so I'm a noob.

I saw this page in documentation which calculates the quadrature using Gauss-Kronrod rule. I was trying to build upon it and calculate a 2-d integral using Gauss-Kronrod rule. My understanding is a 2D integral is basically 2 1D integrals with the Cartesian product. So I was thinking If I can nest two integration function then I'm done. I've started from the aforementioned example and extended it a bit so the code for a 2d integration looks like,

    #include <boost/math/quadrature/gauss_kronrod.hpp>
    #include <iostream>
   
    int main(int argc, char *argv[])
    {
      using namespace boost::math::quadrature;
      auto f1 = [](double t, double s) { return std::exp(-(t*t+s*s+t*s) / 2); };
      double error;
      double Q = gauss_kronrod<double, 15>::integrate(gauss_kronrod<double, 61>::integrate(f1, 0, std::numeric_limits<double>::infinity(), 5), 0, std::numeric_limits<double>::infinity(), 5, 1e-9, &error);
      std::cout << Q << " " << error<<std::endl;
      return 0;
    }
.
I've just added a lambda function (which is the integrand) of two variables and nested 2 1D Gauss-Kronrod integration functions. I'm compiling this with

g++ -Wall -I /path/to/boost_1_73_0  2d-gauss.cpp

This gives error and compilation fails. I don't know what wrong I'm doing. How do you suggest to do such an integral?
Thanks in advance.