Listing 6: A convergence test corresponding to the differential equations shown in Listing 3.

bool approx_equal(const valarray<complex<double> > &x1,
                  const valarray<complex<double> > &x2)
{
    const double rtol2 = 1e-8;
    bool ap_eq(true);
    for (unsigned int loop=0; loop<x1.size(); loop++)
        ap_eq = ap_eq &&
            (norm(x1[loop]-x2[loop])
         < rtol2*norm(x1[loop]+x2[loop]));
    return(ap_eq);
}