A solution with Makie [sorry is not with Plots]
using CairoMakie
f(x,y) = (x + 2y^2) * abs(sin(y) + cos(x))
x = y = 1:0.5:20
z = [f(x,y) for x in x, y in y]
fig = Figure(resolution=(750,300), fontsize = 14)
ax1 = Axis(fig, aspect = 1, xlabel = "x", ylabel = "y")
ax2 = Axis(fig, aspect = 1, xlabel = "x")
ax3 = Axis(fig, aspect = 1, xlabel = "x")
p1 = heatmap!(ax1, x, y, z, colormap = :Spectral)
contour!(ax2, x, y, z, color = :black, levels = 100:1:101) # contourf! also works... but with colormap
heatmap!(ax3, x, y, z, colormap = (:Spectral, 0.3))
contour!(ax3, x, y, z, color = :black, levels = 100:1:101)
cbar = Colorbar(fig, p1, width = 10, hight = Relative(0.8), ticklabelsize = 10)
limits!(ax2, 0,20,0,20)
hideydecorations!(ax2, grid = false)
hideydecorations!(ax3, grid = false)
fig[1,1] = ax1
fig[1,2] = cbar
fig[1,3] = ax2
fig[1,4] = ax3
save("./results/FigMixHeatContour.png", fig, px_per_unit = 2)
fig
