x1 = 0; x2 = 1;
y1 = 0; y2 = 1;
z11 = 0; z12 = 2;
z21 = 3; z22 = 1;
point_n = 21;
x = linspace(x1, x2, point_n);
y = linspace(y1, y2, point_n);
[xx, yy] = meshgrid(x, y);

p1 = (xx-x1)/(x2-x1);
p2 = 1-p1;
q1 = (yy-y1)/(y2-y1);
q2 = 1-q1;
zz = z11*(p2.*q2) + z12*(p2.*q1) + z21*(p1.*q2) + z22*(p1.*q1);

mesh(xx, yy, zz);
line(x1, y1, z11, 'linestyle', '.', 'markersize', 25);
line(x1, y2, z12, 'linestyle', '.', 'markersize', 25);
line(x2, y1, z21, 'linestyle', '.', 'markersize', 25);
line(x2, y2, z22, 'linestyle', '.', 'markersize', 25);
xlabel('X');
ylabel('Y');
zlabel('Z');
set(gca, 'box', 'on');
frot3d on