%------------------------------------------------------------
% Laplacian
%   Given a function, and grid parameters, returns a pointer
%   to the laplacian of that function. The second derivative
%   is computed by the template:
%       1/dx^2 * (f(x+dx) - 2*f(x) + f(x-dx))
%------------------------------------------------------------
% Parameters:
%   f(r): function pointer accepting r = [x y z]
%   gp:   "grid parameters" specifying the
%         integration grid in the following format
%            [x_lower x_upper x_steps;...
%             y_lower y_upper y_steps;...
%             z_lower z_upper z_steps];
%------------------------------------------------------------
function lap_f = laplacian(f,gp)

dx = (gp(1,2) - gp(1,1))/(gp(1,3)-1);
dy = (gp(2,2) - gp(2,1))/(gp(2,3)-1);
dz = (gp(3,2) - gp(3,1))/(gp(3,3)-1);

lap_f = @(r)( 1/dx^2*(   f([r(1)+dx,r(2),r(3)]) - ...
                       2*f([r(1),r(2),r(3)])    + ...
                         f([r(1)-dx,r(2),r(3)]) )       + ...
              1/dy^2*(   f([r(1),r(2)+dy,r(3)]) - ...
                       2*f([r(1),r(2),r(3)])    + ...
                         f([r(1),r(2)-dy,r(3)]))        + ...
              1/dz^2*(   f([r(1),r(2),r(3)+dz]) - ...
                       2*f([r(1),r(2),r(3)])    + ...
                         f([r(1),r(2),r(3)-dz]))               );