TL DR:
$$u_1(x_1) = \cos(2\pi~(\frac{x_1}{L_1}) - \pi) + 1$$
$$u_2(x_2) = \cos(2\pi~(\frac{x_2}{L_2}) - \pi) + 1$$
$$u(x,t) = \exp(-a t) u_1(x_1) u_2(x_1)$$
How to construct it:
Sines and cosines are easily differentiable so they make a good starting point to construct such a solution. We chose a section and offset of the cosine which has a derivative of zero at the domain boundaries (which are 0 and 1 if $L_x=1$):
$$u_1(x_1) = \cos(2\pi~(\frac{x_1}{L_1}) - \pi) + 1$$
(see visualization)
Conveniently the second spatial derivative of $u_1(x_1)$ is:
$$\Delta u_1(x_1) = 4 \pi^2 cos(2 \pi x_1)$$
Now, in higher dimensions, you can construct solutions by multiplying trigonometric functions like this. So you are be able to construct:
$$u(x) = u_1(x_1) u_2(x_2)$$
where u_2 is chosen also as a cosine.
(see 2D visualization)
Finally, the full solution to the heat equation can be guessed by inserting $u(x)$
into:
$$\frac{\partial u}{\partial t} = \Delta u(x,y)$$
When you calculate the l.h.s, you find that it is almost identical to the r.h.s except for a factor. So you may guess that there is an exponential decay of the amplitude, with the variable $t$ as an argument.