

J = 1;
B = 0;
N = 32;

T0 = 1;
T1 = 4;
dT = 0.05;
NT = ceil((T1-T0)/dT)+1;

mPositive = zeros(1,NT);
mNegative = zeros(1,NT);

for n = 1:NT,
    disp(n);
    
    T = T0 + (n-1)*dT;
    beta = 1/T;
    
    [s,Esave,Msave] = ising(beta,J,B,0.99,N);
    mPositive(n) = mean(Msave);
    [s,Esave,Msave] = ising(beta,J,B,0.01,N);
    mNegative(n) = mean(Msave);
    
end
   
plot(T0+[0:NT-1]'*dT,[mPositive;mNegative]','o-');