// From documentation of Asymptote
import graph;
import palette;
import contour;
texpreamble("\usepackage{icomma}");

size(10cm,10cm,IgnoreAspect);

pair a=(0,0);
pair b=(5,10);

real fz(pair z) {
  return z.x*z.y*exp(-z.x);
}
real f(real x, real y) {return fz((x,y));}

int N=200;
int Divs=10;
int divs=2;

defaultpen(1bp);
pen Tickpen=black;
pen tickpen=gray+0.5*linewidth(currentpen);
pen[] Palette=BWRainbow();

scale(false);

bounds range=image(f,Automatic,a,b,N,Palette);
    
xaxis("$x$",BottomTop,LeftTicks(pTick=grey, ptick=lightgrey, extend=true),Above);
yaxis("$y$",LeftRight,RightTicks(pTick=grey, ptick=lightgrey, extend=true),Above);

// Major contours
real[] Cvals;
Cvals=sequence(11)/10 * (range.max-range.min) + range.min;
draw(contour(f,a,b,Cvals,N,operator ..),Tickpen);

// Minor contours
real[] cvals;
real[] sumarr=sequence(1,divs-1)/divs * (range.max-range.min)/Divs;
for (int ival=0; ival < Cvals.length-1; ++ival)
    cvals.append(Cvals[ival]+sumarr);
draw(contour(f,a,b,cvals,N,operator ..),tickpen);

palette("$f(x,y)=xye^{-x}$",range,point(NW)+(0,1),point(NE)+(0,0.25),Top,Palette,
        PaletteTicks(N=Divs,n=divs,Tickpen,tickpen));