function [fx, fy, fz] = gradient_o2(f, dx, dy, dz)
	% calculate second order accurate approximation to gradient of f(x,y,z)
	% using the formula
	% df/dx = (u(i+1,j,k) - u(i-1,j,k))/(2dx)
	
	%% fill in missing arguments
	if nargin < 2
		dx				= 1;
	end
	if nargin < 3
		dy				= dx;
	end
	if nargin < 4
		dz				= dy;
	end
	
	%% generate kernels
	kernel_x			= reshape([1 0 -1] / (2*dx), [3 1 1]);
	kernel_y			= reshape([1 0 -1] / (2*dy), [1 3 1]);
	kernel_z			= reshape([1 0 -1] / (2*dz), [1 1 3]);
	
	%% convolve over given dimension
	fx					= convn(f, kernel_x, 'same');
	fy					= convn(f, kernel_y, 'same');
	fz					= convn(f, kernel_z, 'same');
	fx([1 end],:,:)= nan;
	fy(:,[1 end],:)= nan;
	fz(:,:,[1 end])= nan;
end