

for i=1:6
  beta = 1.5*(i/6)

  n=0:4;
  prob(i,:) = 1./( 1+exp(2*beta*(n-2)) );
end

plot(n,prob)
title('Probability that X_s = 1')
xlabel('Number of neighbors not equal to 1')
ylabel('Probability')
