how can I make this vode work?

clc; 
clear all;
close all;
global phie  e rhocp L dt tf N alpha1 dr drp alpha r
L=0.01; 
dt=0.01; tf=50;
N=21; %M=round(tf/dt);
dr=L/(N-1); drp=2*dr;

alpha=0.5;
r=0.00;
r_out=0.01;
M=21; %nombre de pas entre r_in et r_out
dr=(L-r)/M;
%nr=M+1; %nombre total de noeuds radial
nphi=30;% nombre de pas d'angles
delta_phi=2*pi/nphi; %pas d'angle
t=500; %temps
nt=0.1; % pas dans l'espace temps
dt=t/nt; 
T_in=28;
T_out=28;
err=10^(-3);    Tavant=28; Tarriere=28;Ta=28;

%**********************************  lameff en W/m.K ; I en Ampere ;
%X12=0.02887; c=1495;  K=0.9898;
%rhocp  en J/m.K  X12=4.9054; c=8.4652;  K=0.0082;   
X12=2.887; c=14.95;  K=0.009898;
% inputting the initial values********************      alfa est un coefficient entre 0 et 1 qui dépend de la tortuosité et de la porosité du matériau 
 %X12=4.7419; c=7.1141;  K=0.0082;   %X12,c et K sont les coefficients du modèle de GAB
alpha1=0.5;alpha=0.5; R=8.314; Mv=0.018; rhos=1000; rhoe=382; cps=1240; cpl=4181;% rhos est la masse volumique de bois de tali Ngono 2019 (changer la valeur en fonction de la teneur en eau initiale)
lambda=0.169; rhocp=8.08e5; p=1013.25; pva=1013.25; Lv=2504;% rhocp et lameff correspondent aux valeurs pour le Tali en transversale (Ngono et al 2019)
%********************************** rhos et cps sont respectivement la masse volumiq et la capacité thermik du materiau sec
%********************************** p est la Pression gazeuse totale %lameff la cond thermik effective du matériau
T=ones(N+1);Pvs=ones(N+1);HR=ones(N+1);Da=ones(N+1);pv=ones(N+1);xv=ones(N+1); X12=ones(N+1);c=ones(N+1);K=ones(N+1); 
X=ones(N+1);dHR=zeros(N+1);dPvs=zeros(N+1);DVX=zeros(N+1);DVT=zeros(N+1);lambdae=ones(N+1);gamma=ones(N+1); HRa=ones(N+1); 
T1=ones(M+1,N+1);X1=ones(M+1,N+1);DVX1=zeros(M+1,N+1);DVT1=zeros(M+1,N+1); lambdae1=ones(M+1,N+1); gamma1=ones(M+1,N+1);
%Dm=zeros(N+1);Dm1=zeros(M+1,N+1); Delta=zeros(N+1); Delta1=zeros(M+1,N+1);
 %Ahr=ones(N+1);Bhr=ones(N+1);Chr=ones(N+1); c=ones(N+1);


    re=3949.6159; %Rapport resistance electrique/section sonde (ohm/m2)
    I=0.252;
    phie=(re*I^2)/2;
    rces=585.8002; %rces est la capacité thermique de la sonde
    e=L; % epaisseur du matériau (m)
    lam=lambda; %conductivité isolant (w/m.K)
    roc=rhocp; %capacité isolant (J/m.K);
    %lecture des paramètres
    prompt={'epaisseur materiau (m)','Intensité (A)','conductivité thermique du materiau'};
    def={num2str(e) num2str(I) num2str(lam)};
    titre='Entrée des conditions expérimentales';
    lineNo=[1 60];
    answer=inputdlg(prompt,titre,lineNo,def);

Error using matlab.internal.lang.capability.Capability.require
Support for Java user interfaces is required, which is not available on this platform.

Error in inputdlg (line 55)
Capability.require(Capability.Swing);

    stranswer1=deal(answer{1,:});
    stranswer2=deal(answer{2,:});
    stranswer3=deal(answer{3,:});
    r0=str2double(stranswer1); %str2double plus rapide que str2num mais valide seulement pour les scalaires
    r1=str2double(stranswer2);
    I=str2double(stranswer3);
    %fin lecture des paramètres

for i=2:1:N
T(i)=(Tavant+Tarriere)/2;  Pvs(i)=exp(25.555-(5220/(T(i)+273))); %HR(i)=(HR(1)+HR(N+1))/2;
%teneur en eau initiale calculée par le modèle de GAB
 X(i)=23;
 X12(i)=0.0035*T(i)^2-0.1045*T(i)+4.7088; % X12(i)=-1.3786e-5*T(i)^2+0.0004558*T(i)+0.069;% 
 c(i)=-0.0001*T(i)^2-0.403*T(i)+19.312; % c(i)=0.0064*T(i)^2-0.58*T(i)+22;  %               % le modèle de X12 k et c en fonction de la temperature est tiré de la thèse de SAED RAJI
 K(i)=3.7143e-5*T(i)^2-0.00077*T(i)+0.747; %K(i)=8e-6*T(i)^2-0.0007*T(i)+0.0214;
 Ahr(i)=c(i)-1; Bhr(i)=(X12(i)*c(i))/X(i)-c(i)+2; Chr(i)=1;
%HR(i)=70;
 HR(i)=(1/K(i))*(-Bhr(i)+sqrt(Bhr(i)^2-4*Ahr(i)*Chr(i))/(2*Ahr(i))); 
 pv(i)=HR(i)*Pvs(i);   %pvs la  Pression de vapeur d’eau saturante à la température T
 Da(i)=2.17e-5*((T(i)+273)/273)^1.88; xv(i)=0.622*(pv(i)/(p)); %xv est la Fraction massique de vapeur d’eau 
%X(i)=(X12(i)*c(i)*K(i)*HR(i))/(1-K(i)*HR(i))*(1+c(i)*K(i)*HR(i)-K(i)*HR(i)); %modèle de GAB    
X(i)=(X12(i)*c(i)*K(i)*HR(i))/((c(i)-2)*K(i)*HR(i)-1-(c(i)-1)*K(i)^2*HR(i)^2);
 dHR(i)=(1-K(i)*HR(i))*(1+c(i)*K(i)*HR(i)-K(i)*HR(i))/(X12(i)*c(i)*K(i)+X(i)*K(i)*(1+c(i)*K(i)*HR(i)-K(i)*HR(i))-X(i)*(1-K(i)*HR(i))*(c(i)*K(i)-K(i))); 
 dHR(i)=(1-K(i)*HR(i))*(1+c(i)*K(i)*HR(i)-K(i)*HR(i))/(X12(i)*c(i)*K(i)+X(i)*K(i)*(1+c(i)*K(i)*HR(i)-K(i)*HR(i))-X(i)*(1-K(i)*HR(i))*(c(i)*K(i)-K(i))); 
 dPvs(i)=(5220/(T(i)+273)^2)*Pvs(i);      %dHR(i)=abs(dHR(i));
 DVX(i)=(alpha/rhos*(1-xv(i)))*((Da(i)*Mv)/(R*(T(i)+273)))*Pvs(i)*dHR(i); 
 DVT(i)=(alpha1/rhos*(1-xv(i)))*((Da(i)*Mv)/(R*(T(i)+273)))*HR(i)*dPvs(i);  %rhocp(i)=rhos*(cps+X(i)*cpl);
 lambdae(i)=(lambda+rhos*Lv*DVT(i)); gamma(i)=lambdae(i)/(rhocp);   % lamapp est la conductivité thermique apparente; app la diffusivité apparente du materiau
  %Dm(i)=DVX(i)+DVT(i);Delta(i)=DVT(i)/DVX(i);
  %gamma(i)=lambdae(i)/(rhocp);
%X(i)=0.4061;  %stockage des variables pour t=0
T1(1,i)=T(i)-Ta;
DVX1(1,i)=DVX(i);  lambdae1(1,i)=lambdae(i);
DVT1(1,i)=(DVT(i));  gamma1(1,i)=gamma(i);
X1(1,i)=X(i);      % Dm1(1,i)=Dm(i);
%Delta1(1,i)=Delta(i);
end
% Boundaries conditions 
for i=1
 T(i)=Tavant; X(i)=X(2);
end
T1(1,1)=T(1)-Ta;   X1(1,1)=X(1);
for i=N+1
 T(i)=Tarriere; X(i)=X(N);
end
T1(1,N+1)=T(N+1)-Ta;   X1(1,N+1)=X(N+1);
%loop for time
% rhos=1044;Ta=28;
 time=zeros(M,1);T_new = zeros(N+1,1);  X_new = zeros(N+1,1);%New Temperature calcul.
 DVX_new=zeros(N+1,1); DVT_new=zeros(N+1,1); lambdae_new=zeros(N+1,1); gamma_new=zeros(N+1,1);
%Dm_new=zeros(N+1,1); Delta_new=zeros(N+1,1)

for t=1:1:M  % loop start for time
      %New Temperature calcul.
   for i=2:N    % loop starts for space  
      %ct=((lambda+rhos*Lv*DVT(i+1))/rhocp)*(dt/(dr*dr)); dtt=(gamma(i))*(dt/(dr*dr)); 
      %as=(dt/(dr*dr))*(rhos*Lv/rhocp)*DVX(i+1); bs=(dt/(dr*dr))*(rhos*Lv/rhocp)*DVX(i); 
   %T_new(i)=T(i,t)+ct*T(i+1,t)-(ct+dtt)*T(i,t)+dtt*T(i-1,t)-(rhos*Lv/rhocp)*(as+bs)*X(i,t)+(rhos*Lv/rhocp)*bs*X(i-1,t)+(rhos*Lv/rhocp)*as*X(i+1,t);
   %T_new(i)=T(i,t)-(dtt+ct)*(2+(1/i))*T(i,t)+dtt*T(i-1,t)+(1+(1/i))*ct*T(i+1,t)-(as+bs)*(2+(1/i))*X(i,t)+as*(1+(1/i))*X(i+1,t)+bs*X(i-1,t);
  T_new(i)=T(i,t)+a*((dt/(dr*dr))+(1/r(i))*dt/(dr))*T(i+1,t)-a*(2*(dt/(dr*dr))+(1/r(i))*dt/(dr))*T(i,t)+a*(dt/(dr*dr))*T(i-1,t)+b*((dt/(dr*dr))+(1/r(i))*dt/(dr))*X(i+1,t)-b*(2*(dt/(dr*dr))+(1/r(i))*dt/(dr))*X(i,t)+b*(dt/(dr*dr))*X(i-1,t);

   end
   
  %for the first node  
  T_new(1)=(phie/lambdae(2))*(dr)+(lambdae(2)/lambdae(2))*T_new(2);% +(Lv*rhos/lambdae(1))*(DVX(2)*X(2)-DVX(1)*X(1));
  %for the last node
  T_new(N+1)=Ta;
  T(:,t+1)=T_new;  % temperature update
   %New water content calcul.
  for i=3:N-1
     at=DVX(i+1); bt=DVX(i-1); aw=DVT(i+1)*(dt/(drp*drp)); bw=DVT(i-1)*(dt/(drp*drp)); 
     %X_new(i)=X(i,t)+as*X(i+1,t)-(as+bs)*X(i,t)+bs*X(i-1,t)+at*T(i+1,t)-(at+bt)*T(i,t)+bt*T(i-1,t);
     
     %X_new(i)=X(i,t)+at*(1+(1/i))*X(i+1,t)-(at+bt)*(2+(1/i))*X(i,t)+bt*X(i-1,t)+aw*(1+(1/i))*T(i+1,t)-(as+bw)*(2+(1/i))*T(i,t)+bw*T(i-1,t);
     X_new(i)=X(i,t)+aw*((dt/(dr*dr))+(1/r(i))*dt/(dr))*T(i+1,t)-(aw+bw)*(2*(dt/(dr*dr))+(1/r(i))*dt/(dr))*T(i,t)+bw*(dt/(dr*dr))*T(i-1,t)+at*((dt/(dr*dr))+(1/r(i))*dt/(dr))*X(i+1,t)-(at+bt)*(2*(dt/(dr*dr))+(1/r(i))*dt/(dr))*X(i,t)+bt*(dt/(dr*dr))*X(i-1,t);

  end 
   %for the first node
%    Pvs(1)=exp(25.555-(5220/(T_new(1)+273))); HRa(1)=100*pva/Pvs(1); 
  % X_new(1)=HRa(1)/T_new(1);
 
   X_new(2)=X(2,t)+DVX(3)*(dt/(drp*drp))*(1+(1/3))*X(3,t)-DVX(2)*(dt/(drp*drp))*(2+(1/2))*X(2,t)+DVX(1)*(dt/(drp*drp))*X(1,t)+DVT(3)*(dt/(drp*drp))*(1+(1/3))*T(3,t)-DVT(2)*(dt/(drp*drp))*(2+(1/2))*T(2,t)+DVT(1)*(dt/(drp*drp))*T(1,t);
  
  X_new(1)=X_new(2);
  %for the last node

  X_new(N)=X(N,t)+DVX(N+1)*(dt/(drp*drp))*(1+(1/N+1))*X(N+1,t)-DVX(N)*(dt/(drp*drp))*(2+(1/N))*X(N,t)+DVX(N-1)*(dt/(drp*drp))*X(N-1,t)+DVT(N+1)*(dt/(drp*drp))*(1+(1/N+1))*T(N+1,t)-DVT(N)*(dt/(drp*drp))*(2+(1/N))*T(N,t)+DVT(N-1)*(dt/(drp*drp))*T(N-1,t);
  X_new(N+1)=X_new(N);
   X(:,t+1)=X_new;  %update of water content
   % calculate the new diffusion coefficients
  for i=1:1:N+1
    Pvs(i)=exp(25.555-(5220/(T_new(i)+273)));%HR(i)=0.3;
    X12(i)=-1.3786e-5*T_new(i)^2+0.0004558*T_new(i)+0.069;
    c(i)=0.0064*T_new(i)^2-0.58*T_new(i)+22;
    K(i)=3.7143e-5*T_new(i)^2-0.00077*T_new(i)+0.747;
    Ahr(i)=c(i)-1; Bhr(i)=(X12(i)*c(i))/X_new(i)-c(i)+2; Chr(i)=1;
    HR(i)=(1/K(i))*(-Bhr(i)+sqrt(Bhr(i)^2-4*Ahr(i)*Chr(i))/(2*Ahr(i)));
    Da(i)=2.17e-5*((T_new(i)+273)/273)^1.88; pv(i)=HR(i)*Pvs(i); xv(i)=0.622*(pv(i)/(p)); 
    %X(i)=(X12*c*K*HR(i))/((c-2)*K*HR(i)-1-(c-1)*K^2*HR(i)^2); %On utilise plus ceci car le regime devient dynamique 
    %dHR(i)=(1-K(i)*HR(i))*(1+c(i)*K(i)*HR(i)-K(i)*HR(i))/(X12(i)*c(i)*K(i)+X(i)*K(i)*(1+c(i)*K(i)*HR(i)-K(i)*HR(i))-X(i)*(1-K(i)*HR(i))*(c(i)*K(i)-K(i))); 
    dHR(i)=(1-K(i)*HR(i))*(1+c(i)*K(i)*HR(i)-K(i)*HR(i))/(X12(i)*c(i)*K(i)+X_new(i)*K(i)*(1+c(i)*K(i)*HR(i)-K(i)*HR(i))-X_new(i)*(1-K(i)*HR(i))*(c(i)*K(i)-K(i)));   
    dPvs(i)=(5220/(T_new(i)+273)^2)*Pvs(i);  %dHR(i)=abs(dHR(i));%rhocp(i)=rhos*(cps+X(i)*cpl);
    DVX_new(i)=(alpha/rhos*(1-xv(i)))*((Da(i)*Mv)/(R*(T_new(i)+273)))*Pvs(i)*dHR(i); 
    DVT_new(i)=(alpha1/rhos*(1-xv(i)))*((Da(i)*Mv)/(R*(T_new(i)+273)))*HR(i)*dPvs(i);
    lambdae_new(i)=(lambda+rhos*Lv*DVT_new(i));   gamma_new(i)=lambdae_new(i)/(rhocp);
  end

    %update  new coefficient
    DVX=DVX_new; DVT=DVT_new; 
    lambdae=lambdae_new; gamma=gamma_new;
    %Dm=Dm_new;  Delta=Delta_new;
  for i=1:1:N+1
      T1(t+1,i)=T_new(i)-Ta;        % data stocking
      DVX1(t+1,i)=DVX(i); lambdae1(t+1,i)=lambdae(i);
      DVT1(t+1,i)=(DVT(i)); gamma1(t+1,i)=gamma(i);
      X1(t+1,i)=X_new(i);   %Dm1(t+1,i)=Dm(i);
      %Delta1(t+1,i)=Delta(i);
  end
  %  if abs(v-T)<err
  %     break
  %  end
  %  end    
 end 
X1; T1;

 for t=1:1:M+1
 time(t,1)=t*dt;
 end
 for i = 1:N+1
r(i) =(i)*dr;
end
 

figure('Name','coef de diff non isotherm','NumberTitle','off')
plot(time,DVX1(:,2),'linewidth',5)
%title('Coefficient de diffusion non isotherme'),
xlabel('temps(s)'),ylabel('DVX(m2/s)')
grid on
%subplot(2,2,1)

figure('Name','Coefficient de diffusion isotherme','NumberTitle','off')
plot(time,DVT1(:,2),'linewidth',5)
%title('Coefficient de diffusion isotherme'),
xlabel('time(s)'),ylabel('DVT(m2/s/K)')
grid on
%subplot(2,2,2)

figure('Name','Non isotherm diff coef vs water content','NumberTitle','off')
plot(X1(:,2),DVX1(:,2),'linewidth',5)
%title('Coefficient de diffusion non isotherme sur la face chauffée'),
xlabel('X(%)'),ylabel('DVX(m2/s)')
grid on
%subplot(2,2,3)

figure('Name','isotherm diff coef vs water content','NumberTitle','off')
plot(X1(:,2),DVT1(:,2),'linewidth',5)
%title('Coefficient de diffusion isotherme sur la face chauffée'),
xlabel('X(%)'),ylabel('DVT(m2/s/K)')
grid on
%subplot(2,2,4)

%subplot(1,2,1)
figure('Name','temperature vs lenght','NumberTitle','off')
plot(r,T1(100,:),'-',r,T1(500,:),'-',r,T1(800,:),'-',r,T1(1000,:),'-','linewidth',2.5)
%title('Temperature')%
xlabel('x(m)')
ylabel('T(°C)')
%subplot(1,2,2)
figure('Name','water vs lenght','NumberTitle','off')
plot(r,X1(100,:),'-',r,X1(500,:),'-',r,X1(800,:),'-',r,X1(1000,:),'-','linewidth',2.5)
%title('Teneur en eau')
xlabel('x(m)')
ylabel('X(%)')

%subplot(1,2,1)
figure('Name','temp in wall','NumberTitle','off')
mesh(T)
%title('Temperature dans la matériau')
xlabel('x(m)')
ylabel('temps(s)')
zlabel('T(°C)')
%subplot(1,2,2)
figure('Name','water in wall','NumberTitle','off')
mesh(X)
%title('Teneur en eau dans le matériau')
xlabel('x(m)')
ylabel('temps(s)')
zlabel('X(%)')

Answer 1

[Image Analyst]

In this line

T_new(i)=T(i,t)+a*((dt/(dr*dr))+(1/r(i))*dt/(dr))*T(i+1,t)-a*(2*(dt/(dr*dr))+(1/r(i))*dt/(dr))*T(i,t)+a*(dt/(dr*dr))*T(i-1,t)+b*((dt/(dr*dr))+(1/r(i))*dt/(dr))*X(i+1,t)-b*(2*(dt/(dr*dr))+(1/r(i))*dt/(dr))*X(i,t)+b*(dt/(dr*dr))*X(i-1,t);

you don't yet have "a" defined. But it tries to use "a". You need to assign a value to "a" first before you attempt to use it.