function Q = func1(P00, P01, P10, P11, R000, R001, R010, R011 , R100, R101, R110, R111)

u = [0:0.01:1];
T = [];
S = [];

for u=0:0.01:1
   T = [T; u.^3 u.^2 u 1];
   S = [S; u.^3 u.^2 u 1];
end

T = T';
S = S';

 M = [2 -3 0 1;
      -2 3 0 0;
      1 -2 1 0;
      1 -1 0 0];

 Gx = [P00(1) P01(1) R001(1) R011(1);
       P10(1) P11(1) R101(1) R111(1);
       R000(1) R010(1) 0 0;
       R100(1) R110(1) 0 0];

 Gy = [P00(2) P01(2) R001(2) R011(2);
       P10(2) P11(2) R101(2) R111(2);
       R000(2) R010(2) 0 0;
       R100(2) R110(2) 0 0];

 Gz = [P00(3) P01(3) R001(3) R011(3);
       P10(3) P11(3) R101(3) R111(3);
       R000(3) R010(3) 0 0;
       R100(3) R110(3) 0 0];

%    size(S)
%    size(M)
%    size(Gx)
%    size(T)

Qx = S'*M'*Gx*M*T;
Qy = S'*M'*Gy*M*T;
Qz = S'*M'*Gz*M*T;

surf(Qx, Qy, Qz);

end