Suppose we have a differential equation of the form $\frac{dy}{dx}=f_{p_x,p_y}(y,x)$. After solving this differential equation for particular values of $p_x$ and $p_y$, we obtain $y_{p_x,p_y}(x)$. With this $y_{p_x,p_y}(x)$, we obtain another function $g(y)$, which is obviously dependent on $p_x$, $p_y$ and $x$. So we can express the newly obtained function as $G(p_x,p_y,x)$. Now, we want to integrate $G$ over $p_x$ and $p_y$ to obtain $I(x)=\int_{p_x}\int_{p_y}G(p_x,p_y,x)dp_x\,dp_y$.
However, since the original differential equation does not have an analytic solution, we have to evaluate $I(x)$ numerically. To do this, we have used the following steps:
1.Define a function $f1(p_x,p_y,x)$ that returns $y_{p_x,p_y}(x)$ // using the Runge-Kutta-Fehlberg 78 method.
2.Define a function $f2(x)$ that returns $I(x)$ // using the GNU Scientific Library CQUAD method.
I used C++ for fast computation. It's worth noting that when evaluating the integral $\int_{p_x}\int_{p_y}G(p_x,p_y,x)dp_x\,dp_y$, we have to evaluate $G(p_x,p_y,x)$ for each step in $p_x$ and $p_y$. This means that we have to solve the differential equation for each step of the integration, which can be time-consuming depending on the number of steps taken by $I(x)$. We need to use an appropriate numerical integration method to solve this problem with accuracy in less time. However, since $f1$ is not a general analytic function, and its values depend on the numerical solution obtained using Runge-Kutta-Fehlberg 78, it is difficult to know how the libraries work internally. When CQUAD tries to return value for $f2$, its adaptive algorithm tries to adapt the error within a limit, it may never achieve this limit as $f1$ is not analytic, and the integration might continue indefinitely. This is currently happening, and the code has been running for the past four days without producing a value. Therefore, the question is how to numerically solve this kind of problem and which integration methods to use.