The bivariate equations descibe the surfaces. Are you expecting plots like the following?
Edit: Correcting the dot-product mistakes in the equations.
x1 = linspace(-10, 10, 51);
x2 = linspace( 0, 20, 51);
[X1, X2] = meshgrid(x1, x2);
Y1 = X1.^2 + X1.*X2 - 10;
Y2 = X2 + 3*X1.*(X2.^2) - 57;
figure
surf(X1, X2, Y1), xlabel('x_{1}'), ylabel('x_{2}'), zlabel('y_{1}')
title('Surface y_{1}')
figure
surf(X1, X2, Y2), xlabel('x_{1}'), ylabel('x_{2}'), zlabel('y_{2}')
title('Surface y_{2}')
