What MATLAB tool should I use got coupled partial differential equations?

I'm working on a research project, that is as desiccant wheel system. It include 4 no of coupled partial differential equations variables in single space and time variable. There are 5 no of unknown variables(Ta, Td, Ya, Yd,W), using available auxiliary equations the the number of variables are shorted to 4. There are initial and boundary conditions also to support the solution of partial differential equations. I have inserted the auxiliary equations in MATLAB program and except the unknown variables all other variables will be available from the program. If any other details required than I will provide.

If required then I will provide the range for unknown variables.

Code for other variables is as follows:

% ********** Specification of Desiccant wheel **************
a=0.0018;        % Height of flow passage
b=0.0036;        % Width of flow passage
delta=0.0004;    % Thickness of desiccant layer
epsilon=0.79;    % Porosity in desiccant
sigma=0.7;       % Volume ratio
rho_m=625;       % Matrix material density
k_d=0.175;       % Thermal comductivity of desiccant
c_pm=1030;       % Specific heat of matrix material
c_pd=921;        % Specific heat of Desiccant

% ********** Supporting variable value for governing equations **************
rho_da=1.225;    % Density of dry air
rho_ad=700;      % Density of adsorbent(Silica gel here)
c_pw=4184;       % Specific heat of water
c_pv=1864;       % Specific heat of water at 300 K
c_pda=1005;      % Specific heat of air at constant pressure
P_atm=101325;    % Atmospheric pressure in Pascal 
W_max=0.33;      % Maximum water content in desiccant
R=0.5;           % Separation factor
T_0=256;
P_0=98000;
T=298;

% ********** Operating variables base values **************
u=2;             % Flowing air speed for both sections 

% ********** Auxilary equations *****************
A_csf=a*b/2;                                                                 % Cross-sectional area of flow passage
A_cst=((a+delta)*(b+delta))/2;                                               % Total cross swctional area of one channel
A_csdl=(A_cst-A_csf);                                                        % Cross-sectional of desiccant layer
j=(3+(2*b/a*pi)^2)/(4+(2*b/a*pi)^2);
P_f=b+(2*((b/2)^2+(pi*a/2)^2)^0.5)*j;                                        % Perimeter of flow channel
D_hf=a*(1.0542-0.466*(a/b)-0.1180*(a/b)^2+0.1794*(a/b)^3-0.043*(a/b)^4);     % Hydraulic diameter
Y_d=(0.62188*RH)/((P_atm/P_vs)-RH);                                          % Humidity ratio and relative humidity ratio
P_vs=exp*(23.196-(3816.44/(T_d-46.13)));                                     % Saturation pressure of water vapour
W/W_max==RH/(RH*(1-R));                                                      % Water content in desiccant
q_ads=2800;                                                                  % Heat of adsorption(Kj/kg) - h_fg*(1+0.2843*exp(-10.28*W))
Nu=(Nu_t+Nu_h)/2;                                                            % Nusselt number
Nu_t=1.1791*(1+2.7701*(a/b)-3.1901*(a/b)^2-1.9975*(a/b)^3-0.4966*(a/b)^4);
Nu_h=1.903*(1+0.455*(a/b)+1.2111*(a/b)^2-1.6805*(a/b)^3+0.7724*(a/b)^4-0.1228*(a/b)^5);
rho_d=sigma*rho_ad+(1-sigma)*rho_m;
% From relationship of Sherwood and Nusselt number
Sh=Nu;
D_m=2.303*10^-5*(P_0/P_atm)*(T/T_0)^1.81;
h_m=(Sh*D_m)/D_hf;
h=(Nu*k_d)/D_hf;

[Torsten]

There is no official MATLAB tool available to solve this system of PDEs.

You will have to discretize the equations in space on your own and solve the resulting system of ordinary differential equations by using ODE15S.

Look up "method-of-lines" for details.

If you are willing to accept second-order derivatives for Ya in the first and Ta in the third equation, you could try MATLAB's "pdepe" .