resolution=10;nGrid=2;
h=1.0/(resolution*nGrid);
X = h/2:h:1-h/2;
[x,y,z]=meshgrid(X);

k=1e2;
a=1.0;b=1.0;c = 0.0;
tmp = (c + cos(2*pi.*x) + a.*cos(2*pi.*y) + b.*cos(2*pi.*z));
tmp2 = exp(k.*tmp);
dHda = -(k.*tmp2.*cos(2*pi.*y))./(tmp2+ 1).^2;
dHdb = -(k.*tmp2.*cos(2*pi.*z))./(tmp2+ 1).^2;
dHdc = -(k.*tmp2)./(tmp2+ 1).^2;

disp( '    |  dHda  |  dHdb  |  dHdc  |')
disp(['max |',num2str(max(dHda,[],'all')),'|',num2str(max(dHdb,[],'all')),'|',num2str(max(dHdc,[],'all')),])
disp(['min |',num2str(min(dHda,[],'all')),'|',num2str(min(dHdb,[],'all')),'|',num2str(min(dHdc,[],'all')),])