How to discretize a system of spatial and temporal differential equation for mass and heat transfer whit the Euler method?

I want to make a model for the equations which are found in the attached article, and have some questions regarding them.

Here are the main second order, partial differential equations

And here them after I discretized them.

My 1) question: Is the discretization correct? And if so, than my code correctly describes it? What I mean is for delta t, and delta x, I should use dt and dx (the time and space increments)?

My code:

for j = 2:m-1   % Time Loop
      for i= 2:n-1  % Space Loop
          X(i,j+1) = D*((X(i+1,j)-2*X(i,j)+X(i-1,j))/(dx^2)+2/(r)*(X(i+1,j)-X(i-1,j)/(2*dx)))*dt + X(i,j)
          T(i,j+1) = lambda_p*((T(i+1,j)-2*T(i,j)+T(i-1,j))/(dx^2)+2/(r)*(T(i+1,j)-T(i-1,j)/(2*dx)))*dt/(ro_p*c_p) + T(i,j)
      end  
end

My 2) question is, that some variables like p_sat, contains the result of the Temperature (T).

How to write it correctly? I tried this way, but I'm almost sure, that T(n,m) is not correct, but don't know how to write it.

p_sat = 1e3*10^(7.19621-1730.63/(233.426+(T(n,m)-273.15)));

My 3) question is that boundary conditions look like this

And I'm not sure if I wrote them correctly into my code, like this (the initials condition is that the Mositure contant (X) and the Temperature (T) equals to the initial X and T)?

% Initial and boundary conditions
X  =  zeros(n,m);
T  =  zeros(n,m);
X(1,1) =  X_in; 
T(1,1) =  T_in; 
% Boundary conditions
X(1,:)   =  0;                           
X(n,:)   =  -k*(c_s-c_f)/(D*ro_p);       
T(1,:)   =  0;                                         
T(n,:)   =  -(h*(T(n,m)-T_f) + k*(c_s-c_f)*dH_vap)/lambda_p;    

Here is my full code which of course not runs :(

clear variables
close all
clear all
% 1. Parameters for the space
n = 20;         % number of space steps
Lx = 3.3e-3;    % max space (= radius)
dx = Lx/(n-1);  % space step
x = 0:dx:Lx;    % space increment
% 2. Parameters for the time 
m = 100;        % number of timesteps
tf = 100;         % final time
dt = tf/(m-1);  % time step
t = 0:dt:tf;    % time increment

% Constants
r = 3.3e-3;       % radius of particle (m)
lambda_p = 0.07;  % Thermal conductivity of particle (W/mK)
T_in = 293;       % Initial particle temperature (K)
X_in = 0.082;     % Initial particla moisture content (kg H2O / kg dry matter)
dH_vap = 2.790e6; % Latent heat of vaporization (J/kg)
c_f = 0.01;       % Drying gas water concentration (kg/m3)
T_f = 330;        % Drying air temperature (K)
ro_p = 800;       % Density of particle (kg/m3)
c_p = 1000;       % Heat capacity of particle (J/kgK)   %% ez egy változó amúgy
D = 3e-9;         % Water diffusivity (m2/s)
k = 0.01;         % Mass transfer coefficient (m/s)
h = 400;          % Convective heat transfer  coefficient (W/m2K)
M = 18;           % Molar mass of water (kg/kmol)
R = 8314;         % Universal gas constant (J/kmolK)
a_w = 0.7         % Water activity at the surface (-)


% Initial and boundary conditions
X  =  zeros(n,m);
T  =  zeros(n,m);
X(1,1) =  X_in; 
T(1,1) =  T_in; 

% Variables
p_sat = 1e3*10^(7.19621-1730.63/(233.426+(T(n,m)-273.15)));
c_sat = p_sat*M/(R*T(n,m));
c_s = a_w*c_sat;


% Boundary conditions
X(1,:)   =  0;                           
X(n,:)   =  -k*(c_s-c_f)/(D*ro_p);       
T(1,:)   =  0;                                         
T(n,:)   =  -(h*(T(n,m)-T_f) + k*(c_s-c_f)*dH_vap)/lambda_p;      


% Implementation of the explicit method
for j = 2:m-1   % Time Loop
      for i= 2:n-1  % Space Loop
      
          X(i,j+1) = D*((X(i+1,j)-2*X(i,j)+X(i-1,j))/(dx^2)+2/(r)*(X(i+1,j)-X(i-1,j)/(2*dx)))*dt + X(i,j)
          T(i,j+1) = lambda_p*((T(i+1,j)-2*T(i,j)+T(i-1,j))/(dx^2)+2/(r)*(T(i+1,j)-T(i-1,j)/(2*dx)))*dt/(ro_p*c_p) + T(i,j)
                   
      end  
end

figure (1)
hold all 
for i=1:4:numel(t)
    plot(x,X(:,i),'linewidth',1.5,'DisplayName',sprintf('t = %1.2f',t(i)))
end
a = ylabel('Moisture');
a = xlabel('Radius (m)');
a=title('Moisture change in function of radius');
legend('-DynamicLegend','location','bestoutside');
grid;

figure (2)
[X, T] = meshgrid(x,t);
s2 = mesh(X',T',X);
set(s2,'FaceColor',[0 0 1],'edgecolor',[0 0 0],'LineStyle','--');
a = title('Moisture distribution in the function of Time and Radisus');
a = xlabel('Radius (m)');
a = ylabel('Time (s)');

Sorry for the long question, but I think they can be answered quickly if someone understands programming in MatLab better than I do.

(I tried to use pdepe, but whit it the results were crap)