function [x,y,bcidx] = meshsquare(m1)
% mesh on a square
[x,y] = meshgrid(linspace(0,1,m1));
x = x(:);
y = y(:);
m = m1 ^ 2;
% index to b.c.
bcidx = [1:m1];
for i = 2:m1-1
  bcidx = [bcidx,m1*(i-1)+1,m1*(i-1)+m1];
end
bcidx = [bcidx,m1^2-m1+1:m1^2];
% permutazione random dei nodi (per simulare un brutto ordinamento)
X = [x(:),zeros(m,1)];
Y = [y(:),zeros(m,1)];
X(bcidx,2) = 1;
Y(bcidx,2) = 1;
rand('twister',0) % fissa il generatore random
idx = randperm(m);
X = X(idx,:);
Y = Y(idx,:);
x = X(:,1);
y = Y(:,1);
bcidx = find(X(:,2))'; % row vector
