Fields

[X, Y] = meshgrid(linspace(-2, 2, 140));
Z = X .* exp(-X.^2 - Y.^2);
contourf(X, Y, Z, 16)
colormap coolwarm
hold on

% the gradient, worked out by hand
[Xq, Yq] = meshgrid(linspace(-2, 2, 15));
E = exp(-Xq.^2 - Yq.^2);
quiver(Xq, Yq, (1 - 2*Xq.^2) .* E, -2*Xq.*Yq .* E, 'Color', 'k')
hold off

axis equal
title('f(x, y) = x·e^(−x²−y²) and ∇f')

[X, Y] = meshgrid(linspace(-2.2, 0.8, 360), linspace(-1.25, 1.25, 300));
Zr = zeros(size(X)); Zi = Zr; count = Zr;

for k = 1:60                       % iterate z ← z² + c on the whole grid at once
  Zr2 = Zr.^2 - Zi.^2 + X;
  Zi  = 2*Zr.*Zi + Y;
  Zr  = Zr2;
  count = count + (Zr.^2 + Zi.^2 < 4);
end

imagesc(X(1, :), Y(:, 1), count)
axis xy, axis equal
colormap magma
title('z ↦ z² + c')