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