You will need to use element-wise operators with vector arguments. See Array vs. Matrix Operations for details.
z = @(x1,x2) (1.5-x1+x1.*x2).^2 +(2.25-x1+x1.*(x1.^2))^2+(2.625-x1+x1.*(x2).^3).^2;
x1 = -5:5;
x2 = -5:5;
[X1,X2] = ndgrid(x1,x2);
figure
surf(X1, X2, z(X1,X2))
grid
xlabel('X_1')
ylabel('X_2')
zlabel('z')
syms x1 x2
z = (1.5-x1+x1*x2)^2 +(2.25-x1+x1*(x1^2))^2+(2.625-x1+x1*(x2)^3)^2;
figure
grid
fsurf(z, [-5 5 -5 5])
xlabel('X_1')
ylabel('X_2')
zlabel('z')
EDIT — Corrected typographical errors.
.