


T = 2;
p = 0.1;

J = 1;
B = 0;
N = 100;

beta = 1/T;
[s,Esave,Msave] = ising(beta,J,B,p,N);

subplot(121);
pcolor(s); axis square;
title (['Ising model, \beta=' num2str(beta) ', J=' num2str(J) ', B=' num2str(B)]);

subplot(122);
plot([Esave;Msave]');
set(gca,'xlim',[1 length(Esave)]);
axis square;
title (['Ising model, \beta=' num2str(beta) ', J=' num2str(J) ', B=' num2str(B)]);
legend('Energy/N^2','Spin/N^2','location','east');



