function [fx, fy, fz] = gradient_o4(f, dx, dy, dz)
	% calculate fourth order accurate approximation to gradient of f(x,y,z)
	% using the formula
	% df/dx = (-u(i+2,j,k) + 8u(i+1,j,k) - 8u(i-1,j,k) + u(i-2,j,k))/(12dx)
	% df/dy = (-u(i,j+2,k) + 8u(i,j+1,k) - 8u(i,j-1,k) + u(i,j-2,k))/(12dy)
	% df/dz = (-u(i,j,k+2) + 8u(i,j,k+1) - 8u(i,j,k-1) + u(i,j,k-2))/(12dy)
	
	%% 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 8 0 -8 1] / (12*dx), [5 1 1]);
	kernel_y			= reshape([-1 8 0 -8 1] / (12*dy), [1 5 1]);
	kernel_z			= reshape([-1 8 0 -8 1] / (12*dz), [1 1 5]);
	
	%% convolve over given dimension
	fx					= convn(f, kernel_x, 'same');
	fy					= convn(f, kernel_y, 'same');
	fz					= convn(f, kernel_z, 'same');
	fx([1:2 end-1 end],:,:)		= nan;
	fy(:,[1:2 end-1 end],:)		= nan;
	fz(:,:,[1:2 end-1 end])		= nan;
end