integrals: An amendment to spde()

The iterative version of `spde()` currently assumes that the
coefficient  `a`  of  `Dp`  were invariant in the loop, in which
case the values of the multiplier  `alpha`  are powers of  `a`.
Otherwise the offset  `beta`  of the solution will be incorrect.

I have added two test cases where  `a`  is not invariant.
In the general case,  `alpha`  can be used instead of a power of  `a`
to construct the correct value of  `beta`.
It doesn't hurt to add the test to solve_poly_rde, but it need not be called at all in the case C == 0, since it means that spde() already has produced the solution.