%%%Data%%%
temp=780;
Q=140000;
R=8.3144;
temp_kelvin = temp+273;
d0=0.000041;
ksi=10;

a=-240;
b=912;
temp_p=temp+273;

ga=abs((temp_p-b)/a);

D=d0*exp(-Q/(R * temp_kelvin))*1E10;
dT=0.001;
dx=1;
dx2=dx*dx;


Ddd=D*dT/dx2;

%%%Initialization%%%
T1 = zeros(1, 102);

for i=1:1:10
    T1(i)=0.67;
end
for i=11:1:102
    T1(i)=0.02;
end
T2 = T1;

%%%Computing%%%
for i=1:1:1000000
    for j=2:1:ksi
      T2(j)=(1-2*Ddd)*T1(j)+Ddd*(T1(j-1)+T1(j+1));
    end
    T2(ksi+1)=T2(ksi);
    if T2(ksi)>ga;
        ksi=ksi+1;
        T2(ksi)=0.02;
    end;
    T1 = T2;
end