function [X,tri] = cubesurf(Xb,n)

[A,B]=meshgrid(linspace(0,1,n),linspace(0,1,n));
A=reshape(A,prod(size(A)),1);
B=reshape(B,prod(size(B)),1);
triAB=delaunay(A,B);
A = Xb(1)*ones(size(A)) + (Xb(2)-Xb(1))*A;
B = Xb(1)*ones(size(B)) + (Xb(2)-Xb(1))*B;
X = [ Xb(1)*ones(n^2,1) A B ;
      Xb(2)*ones(n^2,1) A B ;
      A  Xb(1)*ones(n^2,1) B ;
      A  Xb(2)*ones(n^2,1) B ;
      A  B     Xb(1)*ones(n^2,1);
      A  B     Xb(2)*ones(n^2,1)];
  
  tri = [triAB;
         triAB +   n^2*ones(size(triAB));
         triAB + 2*n^2*ones(size(triAB));
         triAB + 3*n^2*ones(size(triAB));
         triAB + 4*n^2*ones(size(triAB));
         triAB + 5*n^2*ones(size(triAB))];
