function retval = plot3D_gnu(f,xmin,xmax,Nx,title1,label,shading_interp)
graphics_toolkit ("gnuplot");
tx = linspace (xmin(1),xmax(1),Nx(1))';
ty = linspace (xmin(2),xmax(2),Nx(2))';
[xx, yy] = meshgrid (tx, ty);
zz=zeros(Nx(1),Nx(2));
for i=1:Nx(1)
  for j=1:Nx(2)
    x1=[xx(i,j) yy(i,j)];
    zz(i,j)=f(x1);
  endfor
endfor
colormap ("default");
[Fx,Fy] = gradient (zz);
surf(xx,yy,zz,Fy+Fy);
if shading_interp!=0
  shading interp;
endif
title(title1);
xlabel(label(1));
ylabel(label(2));
zlabel(label(3));
endfunction
