why the variables in the parfor loop are not sliced or temporary variables in the Monte Carlo simulation?

Dearest Matlab users,

I am dealing for the first time with the parfor command. I have a Monte Carlo simulation code to simulation the materials evolution, it runs a lot of time so I would like to run it faster through multithreading. First I built a matrix contains 40000 sites, for every monte carlo step, 40000 attemps should be gone through.The following is the first my parfor-loop code. when I run the whole programme, error shows Error using ==> parallel funtion in line 598 Error in ==> parallel_function>make_general_channel/channel_general at 894 Index exceeds matrix dimensions.and If I change D2, it says error: The temporary variable D2 in a parfor is uninitialized.or D2 is indexed but not sliced. what should I do to change B2 or D2 to a matrix that can be valid in parfor loop?? Thanks so much!

if true
  N=200;
Q=50;
N_mcs =100;
PT=40;   %total (A+B+C),
PT1=20; % PA<x<=PT1 belongs crystal B, PT1<x<=PT belongs to crystal A
PA=0;   % amorphous C
%%%%%%%%%%%initialization%%%%%%%%%%%%%%%%%%%%%%%%%
rand('state',110)
A1=round(rand(1,N*N)*(PT-1)+1); %A1 or B1 is the matrix of species.
A1=reshape(A1,N,N);
A1=A1';
B1(3:2+N,3:2+N)= A1(1:N,1:N);
B1(3:2+N,1:2)= A1(1:N,N-1:N);
B1(3:2+N,N+3:N+4)= A1(1:N,1:2);
B1(1:2,3:2+N)= A1(N-1:N,1:N);
B1(N+3:N+4,3:2+N)= A1(1:2,1:N);
B1(1:2,1:2)= A1(N-1:N,N-1:N);
B1(1:2,N+3:N+4)= A1(N-1:N,1:2);
B1(N+3:N+4,1:2)= A1(1:2, N-2+1:N);
B1(N+3:N+4,N+3:N+4 )= A1(1:2,1:2);
PercentageC = sum(sum(A1<=PA))/(N*N),%%check the percentage of amorphous material.
PercentageB = sum(sum(PA<A1&A1<=PT1))/(N*N),%%check the percentage of crystal B material.
PercentageA = sum(sum(PT1<A1&A1<=PT))/(N*N),%%check the percentage of crystal A material.
[I,J] = find(A1>PA&A1<=PT1); %find the sites of the crystal B.
NN=length(I);
A2=zeros(N,N);% the sites the amorphous material occupy was assigned to 0.
AA2=round(rand(1,NN)*(Q-1)+1); %A2 or B2 is the matrix of spin number(orientation).
for i=1:NN
A2(I(i),J(i))=AA2(i);
end
[K,L] = find(PT1<A1&A1<=PT); %find the sites of the crystal A.
MM=length(K);
AA3=round(rand(1,MM)*(Q-1)+51); %A2 or B2 is the matrix of spin number(orientation).
for j=1:MM
A2(K(j),L(j))=AA3(j);
end
B2(3:2+N,3:2+N)= A2(1:N,1:N);
B2(3:2+N,1:2)= A2(1:N,N-1:N);
B2(3:2+N,N+3:N+4)= A2(1:N,1:2);
B2(1:2,3:2+N)= A2(N-1:N,1:N);
B2(N+3:N+4,3:2+N)= A2(1:2,1:N);
B2(1:2,1:2)= A2(N-1:N,N-1:N);
B2(1:2,N+3:N+4)= A2(N-1:N,1:2);
B2(N+3:N+4,1:2)= A2(1:2,N-2+1:N);
B2(N+3:N+4,N+3:N+4 )= A2(1:2,1:2);
%rand('state',100)
D2=cell(204,1);
for i=1:204
  D2{i}=B2(i,:)
end
D1=cell(204,1);
for i=1:204
  D1{i}=B1(i,:)
end
m=round(rand(N)*(N-1)+3); 
n=round(rand(N)*(N-1)+3);
s=reshape(m,1,40000);
q=reshape(n,1,40000);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N_mcs=10;
N=200;
for mcs=1:N_mcs %%%How many Monte Carlo snnulahon (MCS) number
parfor jj=1:N*N %%%%%N*N attempts are a MCS
          D2=D2{s(jj)};
          D1=D1{s(jj)};
          bb(jj)=D2(q(jj));
          m1=s(jj);
          n1=q(jj);
          bb=bb(jj);
        if(n1==fix(n1/2)*2)==0 %due to two forms of hexagon. Determine the neighboring-matrix.
           NM1=[m1,n1-1;m1+1,n1;m1,n1+1;m1-1,n1+1;m1-1,n1;m1-1,n1-1];
        else
           NM1=[m1+1,n1-1;m1+1,n1;m1+1,n1+1;m1,n1+1;m1-1,n1;m1,n1-1];
        end
        NB2=[D2{NM1(1,1)} (NM1(1,2)), D2{NM1(2,1)}(NM1(2,2)), D2{NM1(3,1)}(NM1(3,2)),D2{NM1(4,1)}(NM1(4,2)), D2{NM1(5,1)}( NM1(5,2)), D2{NM1(6,1)}( NM1(6,2))];%%orientation number
        NB1=[D1{NM1(1,1)} (NM1(1,2)), D1{NM1(2,1)}(NM1(2,2)), D1{NM1(3,1)}(NM1(3,2)),D1{NM1(4,1)}(NM1(4,2)), D1{NM1(5,1)}( NM1(5,2)), D1{NM1(6,1)}( NM1(6,2))];%%different 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

[Edric Ellis]

I would suggest starting by trying to simplify your code to make a short, self-contained example. http://sscce.org/