where F is a vector/list containing the cell values of f to be calculated, F0 those of the previous time step, V is the volume, deltaT is the time step, and A denotes the "discretized laplacian" (coefficient matrix).

A naive solution of

df/dt = laplacian(laplacian(f))

would be creating an equation of the type

(F - F0)*V/deltaT - A*A*F,

i.e. multiplying the fvMatrix A obtained from fvm::laplacian(f) by itself.

I would be grateful if someone could give me a hint on the implementation or how to solve the coupled problem