No need to dig through the standard library in the hopes of finding a built-in “vectorized” routine for this. A loop (in a function) will be clearer and faster.
function ∇²(f)
∇²f = zero(f) # initialize to zero since we don't touch the boundaries
for y = 2:size(f,2)-1, x = 2:size(f,1)-1 # compute ∇²f in interior of f
∇²f[x,y] = f[x-1,y] + f[x+1,y] + f[x,y-1] + f[x,y+1] - 4f[x,y]
end
return ∇²f
end
Note that I only looped over the “interior” of the array f. You could write separate loops to handle the boundaries. But a much easier solution in this kind of code is to use “ghost cells”: simply define extra array elements around the boundaries and assign them to whatever boundary condition you want (e.g. Dirichlet, periodic, …) before computing the Laplacian.
If you are performing additional calculations on the Laplacian (e.g. explicit timestepping of a heat equation), you can combine the above loop with your other calculations rather than allocating and returning a temporary ∇²f array. There are also fancier tricks using Base.Cartesian and @nloops to write a Laplacian implementation that works for any number of dimensions.
On the other hand, if you are performing sparse solves, e.g. solving Poisson’s equations, then you typically want the Laplacian operator as an explicit sparse matrix. On a Cartesian grid, this can be accomplished neatly with Kronecker products, as shown in this notebook from one of my classes.