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 (1e-15). We also use the algorithm to determine the error curve for a Simpson cubature rule.