@zhk rhop = 10922; rhon = 10922; Cp = 200; Cn = 200; Lan =
1.5*10^(-3); Lap = 1.5*10^(-3); Lbn = 1.5*10^(-3); Lbp = 1.5*10^(-3);

kp = 1.8; kn = 2.2; sigmap = 1/(1.2*10^(-5)); sigman =
1/(10^(-5)); taup = 0.00027; taun = -0.000156; L =
2.325*10^(-3); ha = 10;

Tinf = 298; Th = 300; Tc = 298; A = 2.325*10^(-6); Ic = 0.5;
PDE1 = rhop*Cp*D[Tp[t, x, y, z], t] +
Inactive[Div][-kp Inactive[Grad][Tp[t, x, y, z], {x, y, z}], {x, y,
z}] - 1/sigmap*(Ic/A)^2 + taup*Ic/A*D[Tp[t, x, y, z], x];

PDE2 = rhon*Cn*D[Tn[t, x, y, z], t] +