c = zeros(2000,1000);
c1 = zeros(2000,1000);%change diffusivity
c2 = zeros(2000,1000);%2 diffusivity
x = zeros(2000,1);
cf = zeros(2000,1);
cf1 = zeros(2000,1);
cf2 = zeros(2000,1);
for i = 1:2
    c(i,1) = 1;
    c1(i,1) = 1;
    c2(i,1) = 1;
    x(i,1) = i;
end
for i = 3:2000
    c(i,1) = 0;
    c1(i,1) = 0;
    c2(i,1) = 0;
    x(i,1) = i;
end
for t = 1:49999
    for i = 1: 2
        c(i,t+1) = c(i,t);
    end
    for i = 3:1999
        c(i,t+1) = c(i,t)+ 0.5*(c(i+1,t)+c(i-1,t)-2*c(i,t));
    end
end
for i =1:2000
    cf(i,1) = c(i,40000);
end
% the first diffusivity calculation end (pure M phase)
for t = 1:49999
    for i = 1: 2
        c1(i,t+1) = c1(i,t);
    end
    for i = 3:1999
        c1(i,t+1) = c1(i,t)+ 0.1*(c1(i+1,t)+c1(i-1,t)-2*c1(i,t));
    end
end
for i =1:2000
    cf1(i,1) = c1(i,40000);
end
% the second diffusivity calculation end (pure I phase)
% mix 2 diffusivity
for t = 1:49999
    for i = 1: 2
        c2(i,t+1) = c2(i,t);
    end
    for i = 3:1999
        if c2(i,t) > 0.01 % threshold concentration c0, when c0 ~ 0, the result is close to pure M phase, when c0>>0, the result is close to pure I phase
            c2(i,t+1) = c2(i,t)+ 0.5*(c2(i+1,t)+c2(i-1,t)-2*c2(i,t));
        else
            c2(i,t+1) = c2(i,t)+ 0.1*(c2(i+1,t)+c2(i-1,t)-2*c2(i,t));
        end
    end
end
for i =1:2000
    cf2(i,1) = c2(i,40000);
end
% 2 diffusivity calculation end
plot(x,cf,'LineWidth',1);
hold on 
plot(x,cf1,'LineWidth',1);
hold on
plot(x,cf2,'LineWidth', 1)
%text(1100,0.7,' \leftarrow thicker blue analytical');
%text(1200,0.3,' \leftarrow thinner red numerical');
%text(800,0.5,' D(I) >> D(M)');
%text(800,0.4,' T is away from T(MIT)');

