[I J] = meshgrid(0:6,0:6);
N = I.^2 + I.*J + J.^2;
N = unique(reshape(N, 1, numel(N)));
N = N(N > 7);
N = N(1:15)

