Below are two sample codes for making contour plots of a function of two variables. In the first sample code, the function is expressed as a function of two separate scalar variables x and y.

% plot contours of function f

f = @(x,y) [x^2 + y^2];

% defines grid for x1,x2 each between -4 and 4
x = -4:0.1:4; 
y = -4:0.1:4; 

Nx = length(x);
Ny = length(y);

for i = 1:Nx
  for j = 1:Ny
    F(i,j) = f(x(i), y(j));
  end
end

contourf(x,y,F,[0:20]); % specify contour levels
colorbar
axis equal
axis tight
xlabel('x')
ylabel('y')
title('contours of f(x,y) = x^2 + y^2')

The second script suggests a good initial guess for zeros of the function

i.e. points for which . The script plots contour lines near and . The intersection of these curves are points where both components of are near zero, and so serve as good guesses for a Newton search.

% plot contours of function f

f = @(x) [x(1)^2 + x(2)^2 - 7 ; 1/x(1) - x(2)];

% defines grid for x1,x2 each between 0.1 and 4
x1 = 0.1:0.1:4; 
x2 = 0.1:0.1:4; 

Nx1 = length(x1);
Nx2 = length(x2);

F1 = zeros(Nx1,Nx2);
F2 = zeros(Nx1,Nx2);

% set values of F1, F2 at gridpoints
for i = 1:Nx1
  for j = 1:Nx2

    fx = f([x1(i) ; x2(j)]); % revised version
    F1(i,j) = fx(1);
    F2(i,j) = fx(2);
    
  end
end

figure(1);
hold off
contour(x1,x2,F1,[-0.01 0.01]); % specify contour levels
hold on
contour(x1,x2,F2,[-0.01 0.01]); % specify contour levels
colorbar
axis equal
axis tight
xlabel('x')
ylabel('y')