x = 0:0.01:4;
y = 0:0.01:4;
y2 = 2+y;
z = @(x,y) 16.6698.*(x.*y) - 23.0641.*(x.*y);
plot(x, z(x, y), x, z(x, y2))
Note that this code depends upon x and y being the same shape and orientation. It is not calculating for each combination of x and y, not building a surface.