function [Sij,DEc,dM] = ising2Dpaso(Sij,beta)

%J=1

%extraigo el tamaño de la matriz  Sij de espines en la posicion
%i,j
[Lx, Ly] = size(Sij);
DEc = 0;%accumulated DEc from call 
dM=0; %accumulated dM from call
%recorro cada uno de los sitios para proponer un cambio
for i=1:Lx
  for j=1:Ly
  
  spin_old = Sij(i,j);
  %voy a proponer darlo vuelta
  spin_new = - Sij(i,j); 

  %necesito calcular la energia de la nueva configuracion

  %busco los vecinos del spin ij para ver el cambio de energia, uso la function mod (modulo,
  %"resto") para tener en cuenta las condiciones periodicas.
  iright = mod(i,Lx) + 1; 
  ileft = mod(i+Lx-2,Lx) + 1;

  jup = mod(j,Ly) + 1;
  jdown = mod(j+Ly-2,Ly) + 1;

  DE = 2.0*spin_old*(Sij(iright,j) + Sij(ileft,j) + Sij(i,jup) + Sij(i,jdown));
 


  if ( DE < 0 ) 
        %Acepto
        Sij(i,j) = spin_new;
 	DEc = DEc + DE; 
        dM = dM + (spin_new - spin_old);
  else
    %tiro un numero al azar entre 0 y 1
    p = rand() ; 
    expbetaE = exp(-beta*DE);
    if( p < expbetaE )
    % acepto
     Sij(i,j) = spin_new;
 	DEc = DEc + DE; 
        dM = dM + (spin_new - spin_old);
    end
  end



  end % cierro loop j
end % cierro loop i