Solving moisture diffusion equation

Hi everyone, I am a little bit confused about this,

Now Im stuck in the 'how to set the initial and boundary condition'. and the following figure is what I want to get and the trend should be like the following figure.

and this is my temporary code,

clear; clc;
days=1000;
tcuring=14;
var(1)=0.001;
timeinit=0.001;
counter=2;
while (timeinit<days)
    timeinit=timeinit*1.1;
    var(counter)=timeinit;
    counter=counter+1;
end
timeh=zeros(1,length(var)+1);
timeh(1)=0; %should be changed
for m=1:length(var)
timeh(m+1)=var(m)*24;
end
timeday=timeh/24;
%% Input
c3s=58.6;
c2s=17.2;
c3a=6.97;
c4af=9.9;
Eac3s=30000;    %J/mol
Eac2s=56000;    %J/mol
Eac3a=45000;    %J/mol
Eac4af=32000;   %J/mol
mc=0.4;                           %input mass of cement g/cm3
mw=0.157;                         %input mass of water g/cm3
ms=0.658;                         %input mass of sand g/cm3
mg=1.129;                         %input mass of gravel g/cm3
rh=1;           %relative humidity
%year=100;
days=36500;%*year;      %days
srfarea=330;    %m2/g
refarea=330;    %m2/g
wc=0.5;         %water to cement ratio
To=293.15;      %reference temperature
R=8.314;        %J/mol/K
RH1=linspace(0.005,1,length(timeday));
T=20;           %celcius
Tcal=T+273.15;
rhoeff=2.04;                 %density of HCP at 105C g/cm3
%Nucleation and Growth
Kc3s=1.7;%Maruyama's old paper
Nc3s=0.7;%Lothen
Hc3s=1.4;

Kc2s=0.1;
Nc2s=1;%Lothen
Hc2s=1.55;

Kc3a=1;
Nc3a=0.85;
Hc3a=1.5;

Kc4af=0.37;
Nc4af=0.7;
Hc4af=1.55;

%Diffusion
Kc3sd=1.7;
Kc2sd=0.1;
Kc3ad=0.4;
Kc4afd=0.1;

%Shell formation
Kc3ss=1.7;
Nc3ss=3.3;

Kc2ss=0.1;
Nc2ss=4.1;

Kc3as=2;
Nc3as=2.9;

Kc4afs=0.1;
Nc4afs=3.8;

ccem=0.78;                                          %by lothenbach et al.(2008) and (maruyama et al 2006)
chyd=1.25;                                          %J/gK
cwat=4.1852;
cagg=0.9;                                           %ranging from 0.7-1

Qmax=411.15;%J/g
%% pre-allocating
toalfa=zeros(1,size(var,2)+1);
drc3s=zeros(1,size(var,2)+1);
drc2s=zeros(1,size(var,2)+1);
drc3a=zeros(1,size(var,2)+1);
drc4af=zeros(1,size(var,2)+1);
rc3sn=zeros(1,length(var)+1);
rc3ss=zeros(1,length(var)+1);
rc3sd=zeros(1,length(var)+1);
rc2sn=zeros(1,length(var)+1);
rc2ss=zeros(1,length(var)+1);
rc2sd=zeros(1,length(var)+1);
rc3an=zeros(1,length(var)+1);
rc3as=zeros(1,length(var)+1);
rc3ad=zeros(1,length(var)+1);
rc4afn=zeros(1,length(var)+1);
rc4afs=zeros(1,length(var)+1);
rc4afd=zeros(1,length(var)+1);
%% Parrot-Killoh model
for i=1:length(var) %boundary condition
    if (i==1)
        rc3sn(1)=0;
        rc3ss(1)=0;
        rc3sd(1)=0;
        drc3s(1)=10^-15;
        
        rc2sn(1)=0;
        rc2ss(1)=0;
        rc2sd(1)=0;
        drc2s(1)=10^-15;
        
        rc3an(1)=0;
        rc3as(1)=0;
        rc3ad(1)=0;
        drc3a(1)=10^-15;
        
        rc4afn(1)=0;
        rc4afs(1)=0;
        rc4afd(1)=0;
        drc4af(1)=10^-15;
    else
        rc3sn(1)=0;
        rc3ss(1)=0;
        rc3sd(1)=0;
        drc3s(1)=10^-15;
        
        rc2sn(1)=0;
        rc2ss(1)=0;
        rc2sd(1)=0;
        drc2s(1)=10^-15;
        
        rc3an(1)=0;
        rc3as(1)=0;
        rc3ad(1)=0;
        drc3a(1)=10^-15;
        
        rc4afn(1)=0;
        rc4afs(1)=0;
        rc4afd(1)=0;
        drc4af(1)=10^-15;
    end
end
for j=1:length(var)
    toalfa(j)=(drc3s(j)*c3s+drc2s(j)*c2s+drc3a(j)*c3a+drc4af(j)*c4af)/(c3s+c2s+c3a+c4af);
    %C3S
    if (wc*Hc3s>toalfa(j))
        rc3sn(j+1)=Kc3s/Nc3s*(1-drc3s(j))*(-log(1-drc3s(j)))^(1-Nc3s)*(timeh(j+1)-timeh(j))/24/refarea*srfarea*exp(Eac3s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    else
        rc3sn(j+1)=Kc3s/Nc3s*(1-drc3s(j))*(-log(1-drc3s(j)))^(1-Nc3s)*(timeh(j+1)-timeh(j))/24/refarea*srfarea*(1+3.333*Hc3s*wc-3.333*toalfa(j))^4*exp(Eac3s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    end
    if (wc*Hc3s>toalfa(j))
        rc3ss(j+1)=Kc3ss*(1-drc3s(j))^Nc3ss*(timeh(j+1)-timeh(j))/24*exp(Eac3s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    else
        rc3ss(j+1)=Kc3ss*(1-drc3s(j))^Nc3ss*(timeh(j+1)-timeh(j))/24*(1+3.333*Hc3s*wc-3.333*toalfa(j))^4*exp(Eac3s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    end
    if (wc*Hc3s>toalfa(j))
        rc3sd(j+1)=Kc3sd*(1-drc3s(j))^(2/3)/(1-(1-drc3s(j))^(1/3))*(timeh(j+1)-timeh(j))/24*exp(Eac3s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    else
        rc3sd(j+1)=Kc3sd*(1-drc3s(j))^(2/3)/(1-(1-drc3s(j))^(1/3))*(timeh(j+1)-timeh(j))/24*(1+3.333*Hc3s*wc-3.333*toalfa(j))^4*exp(Eac3s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    end
    %C2S
    if (wc*Hc2s>toalfa(j))
        rc2sn(j+1)=Kc2s/Nc2s*(1-drc2s(j))*(-log(1-drc2s(j)))^(1-Nc2s)*(timeh(j+1)-timeh(j))/24/refarea*srfarea*exp(Eac2s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    else
        rc2sn(j+1)=Kc2s/Nc2s*(1-drc2s(j))*(-log(1-drc2s(j)))^(1-Nc2s)*(timeh(j+1)-timeh(j))/24/refarea*srfarea*(1+3.333*Hc2s*wc-3.333*toalfa(j))^4*exp(Eac2s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    end
    if (wc*Hc2s>toalfa(j))
        rc2ss(j+1)=Kc2ss*(1-drc2s(j))^Nc2ss*(timeh(j+1)-timeh(j))/24*exp(Eac2s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    else
        rc2ss(j+1)=Kc2ss*(1-drc2s(j))^Nc2ss*(timeh(j+1)-timeh(j))/24*(1+3.333*Hc2s*wc-3.333*toalfa(j))^4*exp(Eac2s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    end
    if (wc*Hc2s>toalfa(j))
        rc2sd(j+1)=Kc2sd*(1-drc2s(j))^(2/3)/(1-(1-drc2s(j))^(1/3))*(timeh(j+1)-timeh(j))/24*exp(Eac2s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    else
        rc2sd(j+1)=Kc2sd*(1-drc2s(j))^(2/3)/(1-(1-drc2s(j))^(1/3))*(timeh(j+1)-timeh(j))/24*(1+3.333*Hc2s*wc-3.333*toalfa(j))^4*exp(Eac2s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    end
    %C3A
    if (wc*Hc3a>toalfa(j))
        rc3an(j+1)=Kc3a/Nc3a*(1-drc3a(j))*(-log(1-drc3a(j)))^(1-Nc3a)*(timeh(j+1)-timeh(j))/24/refarea*srfarea*exp(Eac3a/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    else
        rc3an(j+1)=Kc3a/Nc3a*(1-drc3a(j))*(-log(1-drc3a(j)))^(1-Nc3a)*(timeh(j+1)-timeh(j))/24/refarea*srfarea*(1+3.333*Hc3a*wc-3.333*toalfa(j))^4*exp(Eac3a/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    end
    if (wc*Hc3a>toalfa(j))
        rc3as(j+1)=Kc3as*(1-drc3a(j))^Nc3as*(timeh(j+1)-timeh(j))/24*exp(Eac3a/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    else
        rc3as(j+1)=Kc3as*(1-drc3a(j))^Nc3as*(timeh(j+1)-timeh(j))/24*(1+3.333*Hc3a*wc-3.333*toalfa(j))^4*exp(Eac3a/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    end
    if (wc*Hc3a>toalfa(j))
        rc3ad(j+1)=Kc3ad*(1-drc3a(j))^(2/3)/(1-(1-drc3a(j))^(1/3))*(timeh(j+1)-timeh(j))/24*exp(Eac3a/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    else
        rc3ad(j+1)=Kc3ad*(1-drc3a(j))^(2/3)/(1-(1-drc3a(j))^(1/3))*(timeh(j+1)-timeh(j))/24*(1+3.333*Hc3a*wc-3.333*toalfa(j))^4*exp(Eac3a/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    end
    %C4AF
    if (wc*Hc4af>toalfa(j))
        rc4afn(j+1)=Kc4af/Nc4af*(1-drc4af(j))*(-log(1-drc4af(j)))^(1-Nc4af)*(timeh(j+1)-timeh(j))/24/refarea*srfarea*exp(Eac4af/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    else
        rc4afn(j+1)=Kc4af/Nc4af*(1-drc4af(j))*(-log(1-drc4af(j)))^(1-Nc4af)*(timeh(j+1)-timeh(j))/24/refarea*srfarea*(1+3.333*Hc4af*wc-3.333*toalfa(j))^4*exp(Eac4af/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    end
    if (wc*Hc4af>toalfa(j))
        rc4afs(j+1)=Kc4afs*(1-drc4af(j))^Nc4afs*(timeh(j+1)-timeh(j))/24*exp(Eac4af/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    else
        rc4afs(j+1)=Kc4afs*(1-drc4af(j))^Nc4afs*(timeh(j+1)-timeh(j))/24*(1+3.333*Hc4af*wc-3.333*toalfa(j))^4*exp(Eac4af/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    end
    if (wc*Hc4af>toalfa(j))
        rc4afd(j+1)=Kc4afd*(1-drc4af(j))^(2/3)/(1-(1-drc4af(j))^(1/3))*(timeh(j+1)-timeh(j))/24*exp(Eac4af/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    else
        rc4afd(j+1)=Kc4afd*(1-drc4af(j))^(2/3)/(1-(1-drc4af(j))^(1/3))*(timeh(j+1)-timeh(j))/24*(1+3.333*Hc4af*wc-3.333*toalfa(j))^4*exp(Eac4af/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
    end
    %C3S
    if (rc3sn(j+1)<rc3ss(j+1))
        if (rc3sn(j+1)<rc3sd(j+1))
            drc3s(j+1)=drc3s(j)+rc3sn(j+1);
        else
            drc3s(j+1)=drc3s(j)+rc3sd(j+1);
        end
    else
        if (rc3ss(j+1)<rc3sd(j+1))
            drc3s(j+1)=drc3s(j)+rc3ss(j+1);
        else
            drc3s(j+1)=drc3s(j)+rc3sd(j+1);
        end
    end
    %C2S
    if (rc2sn(j+1)<rc2ss(j+1))
        if (rc2sn(j+1)<rc2sd(j+1))
            drc2s(j+1)=drc2s(j)+rc2sn(j+1);
        else
            drc2s(j+1)=drc2s(j)+rc2sd(j+1);
        end
    else
        if (rc2ss(j+1)<rc2sd(j+1))
            drc2s(j+1)=drc2s(j)+rc2ss(j+1);
        else
            drc2s(j+1)=drc2s(j)+rc2sd(j+1);
        end
    end
    %C3A
    if (rc3an(j+1)<rc3as(j+1))
        if (rc3an(j+1)<rc3ad(j+1))
            drc3a(j+1)=drc3a(j)+rc3an(j+1);
        else
            drc3a(j+1)=drc3a(j)+rc3ad(j+1);
        end
    else
        if (rc3as(j+1)<rc3ad(j+1))
            drc3a(j+1)=drc3a(j)+rc3as(j+1);
        else
            drc3a(j+1)=drc3a(j)+rc3ad(j+1);
        end
    end
    %C4AF
    if (rc4afn(j+1)<rc4afs(j+1))
        if (rc4afn(j+1)<rc4afd(j+1))
            drc4af(j+1)=drc4af(j)+rc4afn(j+1);
        else
            drc4af(j+1)=drc4af(j)+rc4afd(j+1);
        end
    else
        if (rc4afs(j+1)<rc4afd(j+1))
            drc4af(j+1)=drc4af(j)+rc4afs(j+1);
        else
            drc4af(j+1)=drc4af(j)+rc4afd(j+1);
        end
    end
%Rate of hydration
rc3s(j+1)=(drc3s(j+1)-drc3s(j))/((timeh(j+1)-timeh(j))/24);
rc2s(j+1)=(drc2s(j+1)-drc2s(j))/((timeh(j+1)-timeh(j))/24);
rc3a(j+1)=(drc3a(j+1)-drc3a(j))/((timeh(j+1)-timeh(j))/24);
rc4af(j+1)=(drc4af(j+1)-drc4af(j))/((timeh(j+1)-timeh(j))/24);
    drc3s(2)=10^-15;
    drc2s(2)=10^-15;
    drc3a(2)=10^-15;
    drc4af(2)=10^-15;
    drc3s(2)=10^-15;
    drc2s(2)=10^-15;
    drc3a(2)=10^-15;
    drc4af(2)=10^-15;
    toalfa(length(var)+1)=(drc3s(length(var)+1)*c3s+drc2s(length(var)+1)*c2s+drc3a(length(var)+1)*c3a+drc4af(length(var)+1)*c4af)/(c3s+c2s+c3a+c4af);
end
for j=1:length(var)
    rtoalfa(j+1)=(toalfa(j+1)-toalfa(j))/((timeh(j+1)-timeh(j))/24);
end
toalfa1=max(toalfa);
%% Heat related properties (assume the volume equals to 1 cm3)
wuc=(1-toalfa)*mc/(mw+mc+ms+mg);               

wf0=mw/(mw+mc+ms+mg);
wf01=mw;
tetahpc=0.23;                                          %based on arosa dabarera(2017)as a minimum water to binder ratio,The value for θhp,c was previously used in many studies for representingchemically bound water of cement (Powers, 1960; Neville, 1995; Lam, Wong, &Poon, 2000; Saengsoy, 2003)
wwhp=tetahpc*mc*toalfa/(mw+mc+ms+mg);                  %chemically bound water content
wwhp1=tetahpc*mc*toalfa;  %g/cm3 blm dalam bentuk rate

tetagelc=0.0126+(0.0026/(-1.009+exp(0.1414*wc)));
wwgel=mc*toalfa*tetagelc/(mw+mc+ms+mg);                %gel water content
wwgel1=mc*toalfa*tetagelc;
wfw=wf0-wwhp-wwgel;                                    %weight of evaporable water (no unit)
wfw1=wf01-wwhp1-wwgel1;
wagg=(ms+mg)/(mw+mc+ms+mg);

whp=1-(wfw+wuc+wagg);

Ccon=wuc*ccem+whp*chyd+wfw*cwat+wagg*cagg;             %J/gK

 
averageQt=toalfa*Qmax;%J/gK
averagedeltat=(averageQt*((c3s+c3a+c2s+c4af)/100)*mc)./((mc+ms+mg+mw)*(Ccon));


vct=0.9;%input(prompt);
if tcuring>5
    if wc<=0.3
        vmbet=(0.068-0.22/tcuring)*(0.85+0.45*0.3)*vct;%2.72*10^-4;
    elseif wc>=0.7
        vmbet=(0.068-0.22/tcuring)*(0.85+0.45*0.7)*vct;
    else
        vmbet=(0.068-0.22/tcuring)*(0.85+0.45*wc)*vct;
    end
else
    if wc<=0.3
        vmbet=(0.068-0.22/5)*(0.85+0.45*0.3)*vct;%2.72*10^-4;
    elseif wc>=0.7
        vmbet=(0.068-0.22/5)*(0.85+0.45*0.7)*vct;
    else
        vmbet=(0.068-0.22/5)*(0.85+0.45*wc)*vct;
    end
end
        
        
Cbet=exp(855/Tcal);
%prompt1='Cement type? type 1=1.1; type 2=1; type 3=1.15; type 4=1.5?';
type=1.1;%input(prompt1);
if tcuring>5
    if wc<=0.3
        n=(2.5+15/tcuring)*(0.33+2.2*0.3)*type;
    elseif wc>=0.7
        n=(2.5+15/tcuring)*(0.33+2.2*0.7)*type;
    else
        n=(2.5+15/tcuring)*(0.33+2.2*wc)*type;
    end
else
    if wc<=0.3
        n=(2.5+15/5)*(0.33+2.2*0.3)*type;
    elseif wc>=0.7
        n=(2.5+15/5)*(0.33+2.2*0.7)*type;
    else
        n=(2.5+15/5)*(0.33+2.2*wc)*type;
    end
end
kbet=(1-1/n)*(Cbet-1)/(Cbet-1);%0.83;
cm=100;                    % Prediction length
h=10;                    % the interval of distance (cm)
M=(cm)*(1/h);             % Cut-off line of x-axis (days)
E=0.0001;
for mm=1:length(RH1)
    if (RH1(mm)<=0.98)
        wrht(mm)=(vmbet*Cbet*kbet*RH1(mm)/((1-kbet*RH1(mm))*(1+(Cbet-1)*kbet*RH1(mm))));
    else
        wrht(mm)=(vmbet*Cbet*kbet*RH1(mm)/((1-kbet*RH1(mm))*(1+(Cbet-1)*kbet*RH1(mm))));
    end
end
for mm=1:length(RH1)
    if (RH1(mm)<=0.98)
        dwrht(mm)=(Cbet*RH1(mm)*kbet^2*vmbet)/((RH1(mm)*kbet*(Cbet - 1) + 1)*(RH1(mm)*kbet - 1)^2) - (Cbet*kbet*vmbet)/((RH1(mm)*kbet*(Cbet - 1) + 1)*(RH1(mm)*kbet - 1)) + (Cbet*RH1(mm)*kbet^2*vmbet*(Cbet - 1))/((RH1(mm)*kbet*(Cbet - 1) + 1)^2*(RH1(mm)*kbet - 1));
    else
        dwrht(mm)=(Cbet*RH1(mm)*kbet^2*vmbet)/((RH1(mm)*kbet*(Cbet - 1) + 1)*(RH1(mm)*kbet - 1)^2) - (Cbet*kbet*vmbet)/((RH1(mm)*kbet*(Cbet - 1) + 1)*(RH1(mm)*kbet - 1)) + (Cbet*RH1(mm)*kbet^2*vmbet*(Cbet - 1))/((RH1(mm)*kbet*(Cbet - 1) + 1)^2*(RH1(mm)*kbet - 1));%50*w0293 - 50*w98ad;
    end
end
ah=1.05-3.8*wc+3.56*(wc)^2;
bh=-14.4+50.4*wc-41.8*(wc)^2;
gamh=31.3-136*wc+162*(wc)^2;
Dh=ah+bh.*(1-2.^(-10.^(gamh*(RH1-1))));
RH=0.766; %relative humidity
RHH=zeros(M,length(var));
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% boundary condition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for J=1:length(var)+1
    RHH(:,1)=1;
    if (J==1)
        RHH(1,J)=0;
    else
        RHH(1,J)=0;
    end
    RHH(M+1,J)=0;                %it's a boundary condition
    
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Crank & Nicolson method %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 r=1/(2*h^2);
 for J=1:length(var)
      for I=2:M
          RHH(I,J+1)=RHH(I,J);
%          %RHH(I+1,J)=RHH(I,J);
      end
     for L=1:500
         LL=L;                   %number of times of calculation
         for I=2:M
             RHH(I,J+1)=(((r*((timeh(J+1)-timeh(J))/(24)))*(1/dwrht(J))*(Dh(J+1)*(RHH(I+1,J+1)-2*RHH(I,J+1)+RHH(I-1,J+1))+Dh(J)*(RHH(I+1,J)-2*RHH(I,J)+RHH(I-1,J))+0))+RHH(I,J));
         end

         S=2*Dh(J)*((timeh(J+1)-timeh(J))/(24))/h^2;
       
     end
 end

depth=0:h:cm;
% figure (1)
% plot(RH1,dwrht)
% grid on
% figure (2)
% plot(RH1,Dh)
% xlabel('RH')
% ylabel('Dh')
% grid on
% 
figure (3)
plot(depth,RHH(:,58),depth,RHH(:,75),depth,RHH(:,123),depth,RHH(:,140),depth,RHH(:,147))
legend('5 hours','24 hours','100 days','500 days','1000 days')
grid on
xlabel('Depth (cm)')
ylabel('RH')

figure (4)
plot(timeday,RHH(2,:),timeday,RHH(3,:),timeday,RHH(4,:),timeday,RHH(5,:),timeday,RHH(6,:),timeday,RHH(7,:))
grid on
legend('10 cm','20 cm','30 cm','40 cm','50 cm')
xlabel('Days')
ylabel('RH')

Hope someone could help me to understand clearly how to perform this equation. Thankyou very much!

[Kevin Isakayoga]

@Bjorn Gustavsson maybe like this? could you help me to set the boundary condition and initial condition? since I dont know where to put that condition. Thank you in advance

[Kevin Isakayoga]

@Image Analyst Help me sensei.

[Bjorn Gustavsson]

You should read up on the Crank-Nicolson method. You still don't form the RHS-matrix and the LHS-matrix and you also don't have the update-step with the inverse of the LHS-matrix.

[Kevin Isakayoga]

I'm sorry, I still dont understand about it :(

I just follow based on example 1D Diffusion

[Bjorn Gustavsson]

Using that example you should end up with matrices AA and BB, and an array d which leads to the C-N-step for solution C:

C(:,i_time+1) = AA\(BB*C(:,i_time) + d)

Where AA and BB are composed of modeling the diffusion and time-variation. Since you have no such step as AA\ or inv(AA) you have not implemented a Crank-Nicolson PDE-integration.

I found it imperative to write out the expressions for the matrix-elements on paper. When I tried to blindly write the code it was a disaster. So I suggest that you write out the C-N matrix-element expressions on paper and then implement it in code. Then try that function on a very small example, on the order of a 10x1 array with constant coefficients for diffusion etc. Then you solve your more difficult problem.

[Kevin Isakayoga]

How about if I use FTCS(Forward Time Centered Space)? @Bjorn Gustavsson

[Bjorn Gustavsson]

Sure you can do that too. Provided you make sure that you keep the time-step small enough to maintain the stability condition. (In my cases I prefer to use C-N to be more robust to that type of instability)

[Kevin Isakayoga]

hmm, in that case, how to check the stability if my increment of time is not constant?

im pretty sure just using ,

stability=Dh*(timeh(J+1)-timeh(J))/(24*h^2);%should be lower than 0.5

since my dt is using while function to get the increment, therefore I'm not sure with my stability calculation.@Bjorn Gustavsson

[Bjorn Gustavsson]

You "simply" have to make sure that you never take longer steps in time than:

If your algorithm suggest taking a longer time-step you have to let the FTCS-scheme do that in sub-steps, and handle everything else accordingly. You don't need to save every step of the solution if that becomes too large, but if dt becomes too long you might get growing oscillations...

[Kevin Isakayoga]

Thank youu for the explanation! Im very grateful to have you. @Bjorn Gustavsson

Oh, I'm sorry to bother you again, I have one last question, do you know how to solve the boundary condition at x=0, t>0 with dH/dx=0 and x=1, t>0 with H=0? I have tried but it seems my thinking wasn't right. I might need from your point of view. Should both of my boundary condition using neumann b.c.?

Thank you very much in advance.