Technically, your code works now. But you get NaN values for some dydt's. The reason might be that you start with a complete 0 vector of initial conditions and thus divide by 0 somwehere. But it's time you start debugging ...
clc
clear
close all
a = 10;
b = 10;
Lp = 2*b + 2*(sqrt(b^2 + (a*pi)^2))*((3 + (2*b/a*pi)^2)/(4 + (2*b/a*pi)^2));
A = 2*a*b;
Pa = 101325; % Pa
Ta = 30; % Celsius
D = (2.256/Pa)*((Ta+273.13)/256)^1.81;
V = 2;
rhoa = 1.15; % kg/m^3
Sha = 2.45;
Ky = rhoa*Sha*D*Lp/4*A;
Nua = 2.45;
ka = 0.0321; % W/(m K)
h = Nua*ka*Lp/4*A;
had = 0.9; % J/kg K (specific heat of adsorption)
fd = 0.914; % kg/m
fm = 12; % kg/m
cd = 0.85; % J/kg K
cw = 4.19; % J/kg K
cm = 0.9; % J/kg K
ca = 1.005; % J/kg K
cv = 2.028; % J/kg K
w10=0.0133; % intital values
Ta10=303;
Td10=303;
RH10=0.5 ;
q10=0.01;
%increment for lenght
Nz=100;
z=linspace(0,1,Nz);
dz=z(2)-z(1);
%time
t=0:3:299;%secs
%inital conditions all zeros
ICi=zeros(1,Nz);
IC=[ICi,ICi,ICi,ICi].';
%6 variables os IC must be 4 c1-c6
%odesolver
[t,y]=ode15s(@(t,y)f(t,y,Lp,A,Pa,Ta,V,rhoa,Ky,ka,h,had,fd,fm,cd,cw,cm,ca,cv,w10,Ta10,Td10,q10,Nz,dz),t,IC);
%re calculation
%define value
q=y(:,1:Nz);
w=y(:,Nz+1:2*Nz);
Td=y(:,2*Nz+1:3*Nz);
Ta=y(:,3*Nz+1:4*Nz);
%bounday cond
q(:,1)=q10; %kg/kgda
w(:,1)=w10; %celcius
Td(:,1)=Td10; %celcius
Ta(:,1)=Ta10;
q(:,end)=w(:,end-1);
w(:,end)=w(:,end-1);
Ta(:,end)=Ta(:,end-1);
Td(:,end)=Td(:,end-1);
%plot
%PDE-ODe
function dydt=f(~,y,Lp,A,Pa,Ta,V,rhoa,Ky,ka,h,had,fd,fm,cd,cw,cm,ca,cv,w10,Ta10,Td10,q10,Nz,dz)
dydt=zeros(length(y),1);
dqdt=zeros(Nz,1);
dwdt=zeros(Nz,1);
dTddt=zeros(Nz,1);
dTadt=zeros(Nz,1);
%define value
q=y(1:Nz);
w=y(Nz+1:2*Nz);
Td=y(2*Nz+1:3*Nz);
Ta=y(3*Nz+1:4*Nz);
%bounday condition
q(1)=q10;%kg/kg
w(1)=w10;%30 celcius
Td(1)=Td10;%30celcius
Ta(1)=Ta10;
q(end)=q(end-1);
w(end)=w(end-1);
Ta(end)=Ta(end-1);
Td(end)=Td(end-1);
%interior
for i=2:Nz-1
pws(i)=exp(23.196-(3816.44/(Td(i)-46.13)));
RH(i)=94.92211.*q(i)^4 - 76.5155.*q(i)^3 + 19.08577.*q(i)^2 + 0.03858.*q(1)+ 0.14278;
weq(i)=0.62188*RH(i)./((Pa./pws(i))-RH(i));
dqdt(i)=-2*Ky*Pa/fd*(weq(i)-w(i)); % one of the sysetm equations but its an ODE
dqdz(i)=(q(i+1)-q(i-1))./2./dz;
dwdz(i)=(w(i+1)-w(i-1))./2./dz;
dTadz(i)=(Ta(i+1)-Ta(i-1))./2./dz;
dTddz(i)=(Td(i+1)-Td(i-1))./2./dz;
dwdt(i)=-V.*dwdz(i)-(fd/2*A*rhoa).*dqdt(i);
dTddt(i)=-((2*h*Lp)/(fd*(cd+q(i)*cw)+fm*cm)).*(Td(i)-Ta(i))-((2*Ky*Lp*had)/(fd*(cd+q(i)*cw)+fm*cm)).*(weq(i)-w(i))-((2*Ky*Lp*cv)/(fd*(cd+q(i)*cw)+fm*cm)).*(weq(i)-w(i));
dTadt(i)=-V.*dTadz(i)-((fd*(cd+q(i).*cw)+fm*cm)/(A*rhoa*(ca+w(i)*cv))).*dTddt(i)+(Ky*Pa*had)/(A*rhoa*(ca+w(i)*cv))*(w(i)-weq(i));
end
dydt=[dqdt;dwdt;dTddt;dTadt];
end
