using Plots
N = 64 # resolution N x 3N
boxsizeX = 0.5
boxsizeY = 1.5
dx = boxsizeX / N
vol = dx^2
xlin = range(0.5*dx, boxsizeX-0.5*dx,length=N)
ylin = range(0.5*dx, boxsizeY-0.5*dx, length=3*N)
ρ=(0*xlin').+(ylin.>0.75).+1
heatmap(ρ)
1 Like