See, for example, this post, which provides a simple N-dimensional implementation of a finite-difference method for a wave equation using a leapfrog explicit timestepping scheme. It would be pretty easy to change it to the diffusion equation, too.
Among other things, it illustrates computing the finite-difference Laplacian in N dimensions cleanly and efficiently using Julia’s Cartesian-indexing code, and the useful technique of “ghost cells” for boundary conditions (which allow you to separate the boundary-condition code from the finite-difference code).
(Whoops, sorry, I see that this is an old thread that was recently revived.)