Evaluating a double integral using Euler's method and Richardson extrapolation

We transform a double integral into a second-order initial value problem, which we solve using Euler's method and Richardson extrapolation. For an example we consider, we achieve accuracy close to machine precision (~10-13). We find that the algorithm is capable of determining the error curve for an arbitrary cubature formula, and we use this feature to determine the error curve for a Simpson cubature rule. We also provide a generalization of the method to the case of nonlinear limits in the outer integral.