% Probablity density
% Note: This is NOT the ground state!
prob = inline('det([ones(size(x)), cos(2*x), sin(x+y), cos(x-y), sin(3*y)])^2','x','y')

% Number of Fermions
Nferm = 5;

% Dimensions
ndim = 2;

% Initial coordinates
x = 2*pi*rand(Nferm,1);
y = 2*pi*rand(Nferm,1);

attempts = 100;
accept = 0;
hopsize = .1;

for i=1:attempts
   xnew = x + hopsize*2*pi*(rand(Nferm,1)-.5);
   ynew = y + hopsize*2*pi*(rand(Nferm,1)-.5);
   if ( rand(1) < (prob(xnew,ynew)/prob(x,y)) )
      x = xnew;
      y = ynew;
      accept = accept + 1;
   end
end
 
% Plot 4 copies of the configuration
plot([x;x+2*pi;x;x+2*pi],[y;y;y+2*pi;y+2*pi],'ko')