I figured this out in Makie.jl and the result is quite stunning (all I had to do was add the arrows, they already had an example of how to plot a surface/wireframe/contour):
using ForwardDiff
using Makie
using StatsMakie
f(x) = -5*x[1]*x[2]*exp(-x[1]^2-x[2]^2)
x = -1:0.1:1.0
y = -1:0.1:1.0
z = [f([i,j]) for i in x, j in y]
# All the x,y coordinates
xy = [(i,j) for i in x for j in y]
# Function to return the gradient at a given point
g = x -> ForwardDiff.gradient(f, x)
# Collect all the gradients and scale down so that arrows aren't too big in the plot
G = 0.05 .* [g([i[1],i[2]]) for i in xy]
# All the x coordinates
xs = [xy[i][1] for i in 1:length(xy)]
# All the y coordinates
ys = [xy[i][2] for i in 1:length(xy)]
# We need u,v for the arrow plot
u = [G[i][1] for i in 1:length(G)]
v = [G[i][2] for i in 1:length(G)]
# Plot in Makie
scene = surface(x,y,z)
xm, ym, zm = minimum(scene_limits(scene))
contour!(scene, x, y, z, levels = 15, linewidth = 2, transformation = (:xy, zm))
arrows!(xs,ys,u,v, arrowsize=0.025, linewidth=0.5, transformation = (:xy, zm))
wireframe!(scene, x, y, z, overdraw = true, transparency = true, color = (:black, 0.1))
center!(scene)
