How to implement boundary condition and diffusion?

Hello Guys,

I was wondering if you could take a look at my code and give me your feedback if it is correct or not. I have some doubts about the results. Also, I don't know how to implement boundary condition and diffusion equations in my code.

The algorithm is attached.

Thank you so much.

% Mu_o & Mu_pr
mu_o = mu_ost * exp(-G/(R*To));
mu_pr = mu_o * G/(R*To^2);

% Initializing Phi & Theta_s

t=0:dt:t_final;
y=dx/2:dx:h-dx/2;
for jj=1:length(y)
    for ii=1:length(t)
        phi=ones(jj,ii);
        theta_s=zeros(jj,ii);
        T=zeros(jj,ii);
    end
end
phi(:,1)=Phi_o;
phi(:,2:length(t))=0;
% Entering Loops
% Defining Temperature

for i=1:length(y)
    for j=1:length(t)
        if t(j)<t1
            T(i,j)=To-(y(i)/h)*delta_T*(1-t(j)/t1);
        else
            T(i,j)=To;
        end
    end
end

% Defining Porosity 

for i=1:length(y)
    for j=1:length(t)
        phi(i,j+1)=dt*(-mu_o*exp(-G/(R*T(i,j))*DPDF))+ phi(i,j);
    end
end
% Defining Theta_f
Kapa=-etha*alpha*B*h^2;

beta=-Kapa*alpha;
theta_f=zeros(length(y)+1,length(t));
Dth_fDt=zeros(length(y)+1,length(t));
for i=2:length(y)
    for j=1:length(t)
        if phi(i,j)> Phi_st
            Dth_fDt=Kapa*(phi(i+1,j)-2*phi(i,j)+phi(i-1,j))/(dx^2)+...
                beta*(theta_f(i+1,j)-2*theta_f(i,j)+theta_f(i-1,j))/(dx^2);
        else
            Dth_fDt=0;
        end
        theta_f(i,j)=theta_f(i,j)+Dth_fDt*dt;
    end
end