clear all

  D = [1 0;0 -1];
  P = [0 1;1 0];
  q = [1/3 1/3];
  J1 = ones(3,1);
  J3 = ones(3,2);
  delta = 2/3;
  B = [0 0;1 0;1 1];
  Bs = (J1*q + (delta/2)*((2*B - J3)*D + J3))*P