  The modular invariant coefficient of the $D^6R^4$ interaction in the low
energy expansion of type~IIB string theory has been conjectured to be a
solution of an inhomogeneous Laplace eigenvalue equation, obtained by
considering the toroidal compactification of two-loop Feynman diagrams of
eleven-dimensional supergravity. In this paper we determine the exact
$SL(2,\mathbb Z)$-invariant solution $f(x+iy)$ to this differential equation
satisfying an appropriate moderate growth condition as $y\to \infty$ (the weak
coupling limit). The solution is presented as a Fourier series with modes
$\widehat{f}_n(y) e^{2\pi i n x}$, where the mode coefficients,
$\widehat{f}_n(y)$ are bilinear in $K$-Bessel functions. Invariance under
$SL(2,\mathbb Z)$ requires these modes to satisfy the nontrivial boundary
condition $ \widehat{f}_n(y) =O(y^{-2})$ for small $y$, which uniquely
determines the solution. The large-$y$ expansion of $f(x+iy)$ contains the
known perturbative (power-behaved) terms, together with precisely-determined
exponentially decreasing contributions that have the form expected of
D-instantons, anti-D-instantons and D-instanton/anti-D-instanton pairs.
