sleep

close all; clear all;

j = [0.5 0; 0 1];
dj = det(j);

n1 = @(e,n) 0.25*(1-e)*(1-n);
n2 = @(e,n) 0.25*(1+e)*(1-n);
n3 = @(e,n) 0.25*(1+e)*(1+n);
n4 = @(e,n) 0.25*(1-e)*(1+n);

f = zeros(8,1);

e = 1/sqrt(3)* [-1 1];
n = 1/sqrt(3)* [1 -1];

x = [1 0 5 0 5 0 1 0];

for i = 1:2
  for j = 1:2
    N = [n1(e(i),n(j)), 0, n2(e(i),n(j)), 0,n3(e(i),n(j)),0,n4(e(i),n(j)),0;
         0,n1(e(i),n(j)), 0, n2(e(i),n(j)), 0,n3(e(i),n(j)),0,n4(e(i),n(j))];

    N_2 = N'*N;
    f = f + N_2*x';
    end
end

disp(f(1:2:end,:))