Metropolis algorithm in MATLAB -->> error with function handles

i have a piece of metropolis algorithm:

mB=5.79*10^(-9); %Bohr magnetone in eV*G^-1
kB=0.86*10^(-4); %Boltzmann in eV*K^-1
%system parameters
L=60; %side square grid
L2=L*L; % total number grid position
Tstep=5; %step in temperature change (K)
Maxstep=10; %max number of steps 
nmcs=5; % cycle numberof Metropolis algorithm

magnet=NaN(1,Maxstep);%store magnetization in "monte carlo images" of sample


%Creation initial point arrangement of magnetic spins
%Outer parameters
H=100000; %Gauss
T=20; % Kelvin

%Energy alteration in spin-reverse
de =@ (i,j) (2*mB*H).*mlat(i,j);
%Metropolis probability
pmetro=@ (i,j) exp(-de(i,j)./(kB*T));

%Creation and display of initial lattice
mlat=2*round(rand(L,L))-1;

mtotal=sum(mlat(:))./L2

% Alteration of system with time
for ii=1:Maxstep
    for imc=1:nmcs
        for i=1:L
            for j=1:L
                 if pmetro(i,j)>=1
                    mlat(i,j)=-mlat(i,j);

                elseif rand<pmetro(i,j)
                    mlat(i,j)=-mlat(i,j);
                end

            end
        end
    end
    magnet(:,ii)=sum(mlat(:))./L2;
    %figure(ii);
    %pcolor(mlat);
   % shading interp;
end


m1=mean(magnet)
error=std(magnet) ./sqrt(numel(magnet))
fprintf('Temperature = %d K',T)

figure(13)
plot(magnet(1,:),'b.')
axis([0 10 0 0.5])
grid on
xlabel('i (Configuration) ')
ylabel('M/(N*mB)')

Now,the problem is in figure(13).The values it gives me are around zero (0.05,0.02..).It supposes to give me values around 0.3.. Generally,the graph its ok,It gives me the right "shape"(it has points) but as i said around zero. I really don't know how to put this post in order to be understood.Maybe i have some mistake in the "magnet"matrix ,i don't know. Anyway,i don't demand from anybody to check it thoroughly ,i am just asking if with a quick look anyone can help.

ΕDIT--->> Also,sometimes when i run the program ,it gives me :

Undefined function or method 'mlat' for input arguments of type 'double'.

Error in ==> @(i,j)(2*mB*H).*mlat(i,j)

Error in ==> @(i,j)exp(-de(i,j)./(kB*T))

Error in ==> metropolis at 39 if pmetro(i,j)>=1

EDIT--->>> I found the "mistake" .In my code i开发者_Python百科n the loops where i have the function "pmetro" i replaced it with the "exp(-(2*mB*H).*mlat(i,j)./(kB*T))" and the program worked just fine!!! Why it didn't work with calling the "pmetro"??How can i overcome this?Is there a problem with function handles in loops?

Blockquote

---

I very strongly suggest that you try writing code without using any function handles until you're really familiar with Matlab.

The line de =@ (i,j) (2*mB*H).*mlat(i,j); is what causes your problems. In Matlab, when you define a function handle that refers to, say, an array, the function handle will use the array as it was at the time of definition. In other words, even though mlat changes inside your loop, mlat(i,j) inside the function de is always the same. In fact, you cannot even run this code unless you have previously defined mlat in the workspace.

You should therefore rewrite the main loop as follows

for iStep = 1:maxStep
   for imc = 1:mcs
      pmetro = $some function of mlat - this can be calculated using the 
                entire array as input
      %# for each element in mlat (and thus pmetro), decide whether 
      %# you have to switch the spin
      switchIdx = pmetro > 1 | pmetro < rand(size(mlat)); 
      mlat(switchIdx) = -mlat(switchIdx);
   end 
   $calculate magnetization$
end

Also, note that there is a command mean to take the average. No need to sum and then divide by the number of elements.