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); alphan = -0.00019; alphap = 0.00015;
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];