An O(h6) method for the interpolation of harmonic functions in rectangular do- mains is described and analysed. The method is based on an earlier cubic spline tech- nique [7], and makes use of recent results concerning the a posteriori correction of interpolatory cubic splines