That is challenging, however it will work with certain restrictions on the matrix dimensions (i.e. they must both be the same sizes) —
f = @(z) 0.07.*z.^2./(0.09+z.^2);
g = @(y,z) 0.003+0.01.*(42./(42+(y-z).^4));
N = 25;
y = linspace(-10, 10, N);
z = linspace(-10, 10, N);
[Y,Z] = ndgrid(y,z);
figure
surf(Y,Z,f(Z))
colormap(turbo)
figure
surf(Y,Z,g(Y,Z))
colormap(turbo)
figure
surf(Y,Z,g(Y,Z)-f(Z))
hold on
contour3(Y,Z,g(Y,Z)-f(Z), [0 0], '-r', 'LineWidth',2) % Plot Contour At 0
hold off
colormap(turbo)
.
