# Three dimensional Laplacian insulated on lateral faces and convectively exposed on transverse faces (updated)

I have the three dimensional Laplacian $$\nabla^2 T(x,y,z)=0$$ representing temperature distribution in a cuboid shaped wall which is exposed to two fluids flowing perpendicular to each other on either of the $$z$$ faces i.e. at $$z=0$$ (ABCD) and $$z=w$$ (EFGH). Rest all the faces are insulated i.e. $$x=0,L$$ and $$y=0,l$$. The following figure depicts the situation.

The boundary conditions on the lateral faces are therefore:

$$-k\frac{\partial T(0,y,z)}{\partial x}=-k\frac{\partial T(L,y,z)}{\partial x}=-k\frac{\partial T(x,0,z)}{\partial y}=-k\frac{\partial T(x,l,z)}{\partial y}=0 \tag 1$$

The bc(s) on the two z-faces are robin type and of the following form:

$$\frac{\partial T(x,y,0)}{\partial z} = p_c\bigg(T(x,y,0)-e^{-b_c y/l}\left[t_{ci} + \frac{b_c}{l}\int_0^y e^{b_c s/l}T(x,s,0)ds\right]\bigg) \tag 2$$

$$\frac{\partial T(x,y,w)}{\partial z} = p_h\bigg(e^{-b_h x/L}\left[t_{hi} + \frac{b_h}{L}\int_0^x e^{b_h s/L}T(x,s,w)ds\right]-T(x,y,w)\bigg) \tag 3$$

$$t_{hi}, t_{ci}, b_h, b_c, p_h, p_c, k$$ are all constants $$>0$$.

I have two questions:

(1) With the insulated conditions mentioned in $$(1)$$ does a solution exist for this system?

(2) Can someone help in solving this analytically ? I tried to solve this using the following approach (separation of variables) but encountered the results which I describe below (in short I attain a trivial solution):

I will include the codes for help:

T[x_, y_, z_] = (C1*E^(γ z) + C2 E^(-γ z))*
Cos[n π x/L]*Cos[m π y/l] (*Preliminary T based on homogeneous Neumann x,y faces *)

tc[x_, y_] =
E^(-bc*y/l)*(tci + (bc/l)*
Integrate[E^(bc*s/l)*T[x, s, 0], {s, 0, y}]);
bc1 = (D[T[x, y, z], z] /. z -> 0) == pc (T[x, y, 0] - tc[x, y]);
ortheq1 =
Integrate[(bc1[[1]] - bc1[[2]])*Cos[n π x/L]*
Cos[m π y/l], {x, 0, L}, {y, 0, l},
Assumptions -> {L > 0, l > 0, bc > 0, pc > 0, tci > 0,
n ∈ Integers && n > 0,
m ∈ Integers && m > 0}] == 0 // Simplify

th[x_, y_] =
E^(-bh*x/L)*(thi + (bh/L)*
Integrate[E^(bh*s/L)*T[s, y, w], {s, 0, x}]);
bc2 = (D[T[x, y, z], z] /. z -> w) == ph (th[x, y] - T[x, y, w]);
ortheq2 =
Integrate[(bc2[[1]] - bc2[[2]])*Cos[n π x/L]*
Cos[m π y/l], {x, 0, L}, {y, 0, l},
Assumptions -> {L > 0, l > 0, bc > 0, pc > 0, tci > 0,
n ∈ Integers && n > 0,
m ∈ Integers && m > 0}] == 0 // Simplify

soln = Solve[{ortheq1, ortheq2}, {C1, C2}];
CC1 = C1 /. soln[[1, 1]];
CC2 = C2 /. soln[[1, 2]];
expression1 := CC1;
c1[n_, m_, L_, l_, bc_, pc_, tci_, bh_, ph_, thi_, w_] :=
Evaluate[expression1];
expression2 := CC2;
c2[n_, m_, L_, l_, bc_, pc_, tci_, bh_, ph_, thi_, w_] :=
Evaluate[expression2];

γ1[n_, m_] := Sqrt[(n π/L)^2 + (m π/l)^2];


I have used Cos[n π x/L]*Cos[m π y/l] considering the homogeneous Neumann condition on the lateral faces i.e. $$x$$ and $$y$$ faces.

Declaring some constants and then carrying out the summation:

m0 = 30; n0 = 30;
L = 0.025; l = 0.025; w = 0.003; bh = 0.433; bc = 0.433; ph = 65.24; \
pc = 65.24;
thi = 120; tci = 30;
Vn = Sum[(c1[n, m, L, l, bc, pc, tci, bh, ph, thi, w]*
E^(γ1[n, m]*z) +
c2[n, m, L, l, bc, pc, tci, bh, ph, thi, w]*
E^(-γ1[n, m]*z))*Cos[n π x/L]*Cos[m π y/l], {n,
1, n0}, {m, 1, m0}];


On executing an plotting at z=0 using Plot3D[Vn /. z -> 0, {x, 0, L}, {y, 0, l}] I get the following:

which is basically 0. On looking further I found that the constants c1, c2 evaluate to 0 for any value of n,m.

More specifically I would like to know if some limiting solution could be developed to circumvent the problem of the constants evaluating to zero

Origins of the b.c.$$2,3$$

Actual bc(s): $$\frac{\partial T(x,y,0)}{\partial z}=p_c (T(x,y,0)-t_c) \tag 4$$ $$\frac{\partial T(x,y,w)}{\partial z}=p_h (t_h-T(x,y,w))\tag 5$$

where $$t_h,t_c$$ are defined in the equation:

$$\frac{\partial t_c}{\partial y}+\frac{b_c}{l}(t_c-T(x,y,0))=0 \tag 6$$ $$\frac{\partial t_h}{\partial x}+\frac{b_h}{L}(t_h-T(x,y,0))=0 \tag 7$$

$$t_h=e^{-b_h x/L}\bigg(t_{hi} + \frac{b_h}{L}\int_0^x e^{b_h s/L}T(x,s,w)ds\bigg) \tag 8$$

$$t_c=e^{-b_c y/l}\bigg(t_{ci} + \frac{b_c}{l}\int_0^y e^{b_c s/l}T(x,s,0)ds\bigg) \tag 9$$

It is known that $$t_h(x=0)=t_{hi}$$ and $$t_c(y=0)=t_{ci}$$. I had solved $$6,7$$ using the method of integrating factors and used the given conditions to reach $$8,9$$ which were then substituted into the original b.c.(s) $$4,5$$ to reach $$2,3$$.

Non-dimensional formulation The non-dimensional version of the problem can be written as:

(In this section $$x,y,z$$ are non-dimensional; $$x=x'/L,y=y'/l,z=z'/w, \theta=\frac{t-t_{ci}}{t_{hi}-t_{ci}}$$)

Also, $$\beta_h=h_h (lL)/C_h, \beta_c=h_c (lL)/C_c$$ (However, this information might not be needed)

$$\lambda_h \frac{\partial^2 \theta_w}{\partial x^2}+\lambda_c \frac{\partial^2 \theta_w}{\partial y^2}+\lambda_z \frac{\partial^2 \theta_w}{\partial z^2}=0 \tag A$$

In $$(A)$$ $$\lambda_h=1/L^2, \lambda_c=1/l^2, \lambda_z=1/w^2$$

$$\frac{\partial \theta_h}{\partial x}+\beta_h (\theta_h-\theta_w) = 0 \tag B$$

$$\frac{\partial \theta_c}{\partial y} + \beta_c (\theta_c-\theta_w) = 0 \tag C$$

The z-boundary condition then becomes:

$$\frac{\partial \theta_w(x,y,0)}{\partial z}=r_c (\theta_w(x,y,0)-\theta_c) \tag D$$ $$\frac{\partial \theta_w(x,y,w)}{\partial z}=r_h (\theta_h-\theta_w(x,y,w))\tag E$$

$$\theta_h(0,y)=1, \theta_c(x,0)=0$$

Here $$r_h,r_c$$ are non-dimensional quantities ($$r_c=\frac{h_c w}{k}, r_h=\frac{h_h w}{k}$$).

• The thing I notice is if you look at the simplified result of ortheq1, the only way it can be true is if C1 + C2 = 0 and C1 - C2 = 0 which of course means that both must be zero. It doesn't matter what happens after that, so you might examine the bc going into ortheq1 Jul 24, 2020 at 5:41
• @BillWatts appreciate the observation. I too had a similar doubt because while solving the system by hand I encountered a 0=0 situation while applying the orthogonality. This happens because of the Cos[n π x/L]*Cos[m π y/l] functions in the preliminary T distribution. So when we multiply with Cos[k π x/L]*Cos[j π y/l] on both sides and integrate the integrals vanish. This makes me doubt my preliminary T distribution. Can you comment on whether the form T[x_, y_, z_] = (C1*E^(γ z) + C2 E^(-γ z))*Cos[n π x/L]*Cos[m π y/l] is a correct assumption ? I hope I could make myself clear. Jul 24, 2020 at 6:57
• It certainly seems that the system should have a unique solution because the net flux across all the boundaries is zero (homogeneous Neumann on x,y faces and energy transferred from one to the other fluid is equal through the z faces). Additionally, after reading your comment I tried the problem with swapped signs in the RHS of bc1 but came to the same result. Jul 24, 2020 at 6:59
• @IndrasisMitra Why you are talking about analytical solution while try to get numerical solution in the form of series? Jul 24, 2020 at 16:18
• @AlexTrounev Sir, probably my choice of words have been wrong. I was aiming for a solution using the method of separation of variables whose results are in the form of series. In other words, I am aiming to evaluate the three temperature fields $T,t_h,t_c$ in the form of a series solution. Jul 24, 2020 at 16:30

Problem Statement

For notational simplicity, use the non-dimensional formulation described near the end of the question. (Doing so also facilitates comparison with results of a 2D approximation solved earlier.) The PDE is given by

λh D[θw[x, y, z], x, x] + λc D[θw[x, y, z], y, y] + λz D[θw[x, y, z], z, z] == 0


over the domain {x, 0, 1}, {y, 0, 1}, {z, 0, 1}. Normal derivatives vanish on the x and y boundaries. Conditions on the z boundaries are given by

(D[θw[x, y, z], z] + rh (θw[x, y, z] - θwh[x, y]) == 0) /. z -> 1
(D[θw[x, y, z], z] - rh (θw[x, y, z] - θwc[x, y]) == 0) /. z -> 0


with θwc and θwh specified by

D[θwh[x, y], x] + bh (θw[x, y, 1] - θwh[x, y]) == 0
θwh[0, y] == 1
D[θwc[x, y], y] + bc (θw[x, y, 0] - θwh[x, y]) == 0
θwc[x, 0] = 0


Although the solution of the PDE itself can be expressed as a sum of trigonometric functions, the z boundary conditions couple what otherwise would be separable eigenfunctions.

Coupling Coefficients

The coupling coefficients in question are given by

DSolveValue[{D[θc[y], y] + b (θc[y] - 1) == 0, θc[0] == 0}, θc[y], y] // Simplify
a00 = Simplify[Integrate[% , {y, 0, 1}]]
an0 = Simplify[Integrate[%% 2 Cos[n π y], {y, 0, 1}], Assumptions -> n ∈ Integers]
(* 1 - E^(-b y) *)
(* (-1 + b + E^-b)/b *)
(* -((2 b E^-b ((-1)^(1 + n) + E^b))/(b^2 + n^2 π^2)) *)

DSolveValue[{D[θc[y], y] + b (θc[y] - Cos[m Pi y]) == 0, θc[0] == 0}, θc[y], y] // Simplify
a0m = Simplify[Integrate[%, {y, 0, 1}], Assumptions -> m ∈ Integers]
amm = Simplify[Integrate[%% 2 Cos[m π y], {y, 0, 1}], Assumptions -> m ∈ Integers]
anm = FullSimplify[Integrate[%%% 2 Cos[n π y], {y, 0, 1}], Assumptions -> (m | n) ∈ Integers]
(* (b (-b E^(-b y) + b Cos[m π y] + m π Sin[m π y]))/(b^2 + m^2 π^2) *)
(* (b ((-1)^(1 + m) + E^-b))/(b^2 + m^2) *) *)
(* (b^2 E^-b (b^2 E^b + 2 b ((-1)^m - E^b) + E^b m^2 π^2))/(b^2 + m^2 π^2)^2 *)
(* (E^-b (2 (-1)^n b^3 (m - n) (m + n) + 2 b E^b (n^2 (b^2 + m^2 π^2) + (-1)^(1 + m + n)
m^2 (b^2 + n^2 π^2))))/((m - n) (m + n) (b^2 + m^2 π^2) (b^2 + n^2 π^2)) *)

a[nn_?IntegerQ, mm_?IntegerQ] := Which[nn == 0 && mm == 0, a00, mm == 0, an0, nn == 0, a0m,
nn == mm, amm, True, anm] /. {n -> nn, m -> mm}


General Solution

Express the solution as a sum of eigenfunctions in the absence of the z boundary conditions.

λh D[θw[x, y, z], x, x] + λc D[θw[x, y, z], y, y] + λz D[θw[x, y, z], z, z];
Simplify[(% /. θw -> Function[{x, y, z}, Cos[nh Pi x] Cos[nc Pi y] θwz[z]])/
(Cos[nh Pi x] Cos[nc Pi y])] /. π^2 (nc^2 λc + nh^2 λh) -> k[nh, nc]^2 λz
Flatten@DSolveValue[% == 0, θwz[z], z] /. {C[1] -> c1[nh, nc], C[2] -> c2[nh, nc]}
(* -λz k[nh, nc]^2 θwz[z] + λz (θwz''[z] *)
(* E^(z k[nh, nc]) c1[nh, nc] + E^(-z k[nh, nc]) c2[nh, nc] *)


as expected. Note that k[nh, nc] = Sqrt[π^2 (nc^2 λc + nh^2 λh)/λz] has been introduced for notational simplicity. Although the result above for θwz is correct, rearranging the terms a bit is helpful in what follows.

sz = c2[nh, nc] Sinh[k[nh, nc] z]/Cosh[k[nh, nc]] +
c1[nh, nc] Sinh[k[nh, nc] (1 - z)]/Sinh[k[nh, nc]];


because sz /. z -> 0, needed in the z = 0 boundary condition, reduces to c1[nh, nc]. Next, use that boundary condition to eliminate c2[nh, nc].

sθc[nh_?IntegerQ, nc_?IntegerQ] := Sum[a[nc, m] c1[nh, m], {m, 0, maxc}] /. b -> bc;
(D[sz, z] == rc (sz - sθc[nh, nc])) /. z -> 0;
Solve[%, c2[nh, nc]] // Flatten // Apart;
sz1 = Simplify[sz /. %] // Apart
(* (c1[nh, nc] (Cosh[z k[nh, nc]] k[nh, nc] + rc Sinh[z k[nh, nc]]))/k[nh, nc]
- (rc Sinh[z k[nh, nc]] sθc[nh, nc])/k[nh, nc] *)


Finally, use the z = 1 boundary condition to produce a matrix equation for c[nh, nc].

szθh[nh_?IntegerQ, nc_?IntegerQ] := Evaluate[sz1 /. z -> 1]
sθh[nh_?IntegerQ, nc_?IntegerQ] := Evaluate[Sum[a[nh, m] szθh[m, nc], {m, 0, maxh}]]
eq = Simplify[(D[sz1, z] + rh (sz1 - sθh[nh, nc])) /. z -> 1] -
rh (DiscreteDelta[nh] - a[nh, 0]) DiscreteDelta[nc]
(* -rh DiscreteDelta[nc] (-a[nh, 0] + DiscreteDelta[nh]) + (1/k[nh, nc])
(c1[nh, nc] ((rc + rh) Cosh[k[nh, nc]] k[nh, nc] + rc rh Sinh[k[nh, nc]] +
k[nh, nc]^2 Sinh[k[nh, nc]]) - rc rh Sinh[k[nh, nc]] sθc[nh, nc] -
k[nh, nc] (rc Cosh[k[nh, nc]] sθc[nh, nc] + rh sθh[nh, nc])) *)


The source term rh (DiscreteDelta[nh] - a[nh, 0]) DiscreteDelta[nc] arises from θwh[0, y] == 1 instead of equaling zero.

Specific solution for Parameters Set to Unity.

maxh = 3; maxc = 3; λh = 1; λc = 1; λz = 1; bh = 1; bc = 1; rh = 1; rc = 1;
ks = Flatten@Table[k[nh, nc] -> Sqrt[π^2 (nc^2 λc + nh^2 λh)/λz],
{nh, 0, maxh}, {nc, 0, maxc}]
eql = N[Collect[Flatten@Table[eq /. Sinh[k[0, 0]] -> k[0, 0], {nh, 0, maxh}, {nc, 0, maxc}]
/. b -> bh, _c1, Simplify] /. ks] /. c1[z1_, z2_] :> Rationalize[c1[z1, z2]];
(* {k[0, 0] -> 0, k[0, 1] -> π, k[0, 2] -> 2 π, k[0, 3] -> 3 π, k[1, 0] -> π,
k[1, 1] -> Sqrt[2] π, k[1, 2] -> Sqrt[5] π, k[1, 3] -> Sqrt[10] π, k[2, 0] -> 2 π,
k[2, 1] -> Sqrt[5] π, k[2, 2] -> 2 Sqrt[2] π, k[2, 3] -> Sqrt[13] π,
k[3, 0] -> 3 π, k[3, 1] -> Sqrt[10] π, k[3, 2] -> Sqrt[13] π, k[3, 3] -> 3 Sqrt[2] π} *)


eql the numericized version of eq is too long to reproduce here. And, trying to solve eq itself is far too slow. Next, compute the c1 and from them the solution.

Union@Cases[eql, _c1, Infinity];
coef = NSolve[Thread[eql == 0], %] // Flatten
sol = Total@Simplify[Flatten@Table[sz1 Cos[nh Pi x] Cos[nc Pi y] /.
Sinh[z k[0, 0]] -> z k[0, 0], {nh, 0, maxh}, {nc, 0, maxc}], Trig -> False] /. ks /. %;
(* {c1[0, 0] -> 0.3788, c1[0, 1] -> -0.0234913, c1[0, 2] -> -0.00123552,
c1[0, 3] -> -0.00109202, c1[1, 0] -> 0.00168554, c1[1, 1] -> -0.0000775391,
c1[1, 2] -> -5.40917*10^-6, c1[1, 3] -> -4.63996*10^-6, c1[2, 0] -> 4.19045*10^-6,
c1[2, 1] -> -1.24251*10^-7, c1[2, 2] -> -1.17696*10^-8, c1[2, 3] -> -1.02576*10^-8,
c1[3, 0] -> 1.65131*10^-7, c1[3, 1] -> -3.41814*10^-9, c1[3, 2] -> 3.86348*10^-10,
c1[3, 3] -> -3.48432*10^-10} *)


Here are several plots of the solution, beginning with a 3D contour plot.

ContourPlot3D[sol, {x, 0, 1}, {y, 0, 1}, {z, 0, 1}, Contours -> {.4, .5, .6},
ContourStyle -> Opacity[0.75], PlotLegends -> Placed[Automatic, {.9, .9}],
ImageSize -> Large, AxesLabel -> {x, y, z}, LabelStyle -> {15, Bold, Black}]


Next are slices through the solution at the ends and mid-point in z. The second, at z = 1/2, is similar to the seventh plot in the 2D thin slab approximation, even though the calculation here is for a cube.

Plot3D[sol /. z -> 0, {x, 0, 1}, {y, 0, 1}, ImageSize -> Large,
AxesLabel -> {x, y, "θw(z=0)"}, LabelStyle -> {15, Bold, Black}]


 Plot3D[sol /. z -> 1/2, {x, 0, 1}, {y, 0, 1}, ImageSize -> Large,
AxesLabel -> {x, y, "θw(z=0)"}, LabelStyle -> {15, Bold, Black}]


 Plot3D[sol /. z -> 1, {x, 0, 1}, {y, 0, 1}, ImageSize -> Large,
AxesLabel -> {x, y, "θw(z=0)"}, LabelStyle -> {15, Bold, Black}]


Finally, here are θwc and θwh, each computed in two distinct ways, by the expansion given above and by direct integration using the expansion of θw. (They differ in that the latter does not employ the a matrix.) The two methods agree very well except at the edges in y and x, respectively, where convergence of the cosine series is nonuniform. Increasing the number of modes reduces this modest disagreement but does not change sol by more than 10^-4.

Simplify[(sol - D[sol, z]/rc) /. z -> 0, Trig -> False];
DSolveValue[{θwc'[y] + bc (θwc[y] - sol /. z -> 0) == 0, θwc[0] == 0},
θwc[y], {y, 0, 1}] // Chop;
Plot3D[{%, %%}, {x, 0, 1}, {y, 0, 1}, ImageSize -> Large,
AxesLabel -> {x, y, θwc}, LabelStyle -> {15, Bold, Black}]


Simplify[(sol + D[sol, z]/rh) /. z -> 1, Trig -> False];
DSolveValue[{θwh'[x] + bh (θwh[x] - sol /. z -> 1) == 0, θwh[0] == 1},
θwh[x], {x, 0, 1}] // Chop;
Plot3D[{%, %%}, {x, 0, 1}, {y, 0, 1}, ImageSize -> Large,
AxesLabel -> {x, y, θwh}, LabelStyle -> {15, Bold, Black}]


The computation shown here required only a few minutes on my PC.

• @IndrasisMitra I used the three λ merely to rescale the dimensions, so that x, y, and z would vary from zero to one. Note that doing so also required rescaling rc, rh, bc, and bh. By the way, for some parameters it should be possible to obtain approximate symbolic solutions. Dec 29, 2020 at 13:23
• @IndrasisMitra Although I do not have time now to run various other cases, I did correct a scaling error in the last two plots. Dec 31, 2020 at 15:54
• @IndrasisMitra The source term is Exp[-b x]. A Fourier Cosine decomposition of this expression yields the expression rh (DiscreteDelta[nh] - a[nh, 0]) DiscreteDelta[nc]. Hope this helps. By the way, is is good practice to delete comments that no longer contribute to the discussion. Best wishes. Jan 1, 2021 at 22:57
• @IndrasisMitra Compare the two expressions term by term, and you will find that the coefficients are identical. Jan 17, 2021 at 15:55
• @IndrasisMitra Yes, these are orthogonality relations. Check any text on discrete Fourier transforms. The last two plots above are measures of their accuracy. The series converge, but not uniformly. The are less accurate at the end points. Jan 28, 2021 at 16:06

We can solve this problem with using method explained in my answer here and in my paper attached to this post. We solve in the unit cube system of equations

eq1 = λh D[θw[x, y, z], x,
x] + λc D[θw[x, y, z], y,
y] + λz D[θw[x, y, z], z, z] ==
0; bc1 = {(D[θw[x, y, z], z] +
rh (θw[x, y, z] - θwh[x, y]) == 0) /. z -> 1,
(D[θw[x, y, z], z] -
rc (θw[x, y, z] - θwc[x, y]) == 0) /.
z -> 0}; eq2 =
D[θwh[x, y], x] +
bh (θw[x, y, 1] - θwh[x, y]) == 0;
bc2 = θwh[0, y] == 1;
eq3 = D[θwc[x, y], y] -
bc (θw[x, y, 0] - θwc[x, y]) == 0;
bc3 = θwc[x, 0] == 0;


First we generate base functions and solution to the problem as follows

E[m_, t_] := Cos[m t] Exp[-m t]
nn = 5;
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; ycol =
xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2,
nn + 1}]; zcol = ycol; Psijk =
Table[UE[n, t1], {n, 0, nn - 1}]; Int1 = Integrate[Psijk, t1];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y; int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y; M = nn;
M = nn; U1 = Array[a1, {M, M, M}]; U2 = Array[a2, {M, M, M}]; U3 =
Array[a3, {M, M, M}]; B1 = Array[b1, {M, M}]; B2 =
Array[b2, {M, M}]; B3 = Array[b3, {M, M}]; G1 =
Array[g1, {M, M}]; G2 = Array[g2, {M, M}]; G3 =
Array[g3, {M, M}]; G4 = Array[g4, {M, M}]; G5 = Array[g5, {M, M}];
H1 = Array[h1, {M}]; H2 = Array[h2, {M}];

thx[x_, y_] := (Psi[x] . G5 . Psi[y]);
tcy[x_, y_] := (Psi[x] . G4 . Psi[y]);
th[x_, y_] := (int1[x] . G5 . Psi[y]) + H2 . Psi[y];
tc[x_, y_] := (Psi[x] . G4 . int1[y]) + H1 . Psi[x];

u1[x_, y_, z_] := (int2[x] . U1 . Psi[y]) . Psi[z] +
x Psi[y] . G1 . Psi[z] + Psi[y] . B1 . Psi[z];
u2[x_, y_, z_] := (Psi[x] . U2 . int2[y]) . Psi[z] +
y Psi[x] . G2 . Psi[z] + Psi[x] . B2 . Psi[z];
u3[x_, y_, z_] := (Psi[x] . U3 . Psi[y]) . int2[z] +
z Psi[x] . G3 . Psi[y] + Psi[x] . B3 . Psi[y];
uz[x_, y_, z_] := (Psi[x] . U3 . Psi[y]) . int1[z] +
Psi[x] . G3 . Psi[y];
uy[x_, y_, z_] := (Psi[x] . U2 . int1[y]) . Psi[z] +
Psi[x] . G2 . Psi[z];
ux[x_, y_, z_] := (int1[x] . U1 . Psi[y]) . Psi[z] +
Psi[y] . G1 . Psi[z];
uxx[x_, y_, z_] := (Psi[x] . U1 . Psi[y]) . Psi[z];
uyy[x_, y_, z_] := (Psi[x] . U2 . Psi[y]) . Psi[z];
uzz[x_, y_, z_] := (Psi[x] . U3 . Psi[y]) . Psi[z];


Parameters of the model, equations on the grid and variables definition

(*Another set of parameters can be \
bh=bc=2.065,rh=rc=0.861,\[Lambda]x=\[Lambda]y=0.0118,\[Lambda]z=0.\
8162.These parameters correspond to a miniaturized steel (k=16W/mK) \
wall where L=l=25 mm,w=3 mm with water (cp=4178 J/kgK) flowing on \
either side with a mass flow rate of 0.9775 gm/sec.The heat transfer \
coefficient (h) is set to 4590 W/sq.m K.*)
bh = bc = 2.065; rh =
rc = 0.861; λh = λc = 0.0118; λz = 0.8162;

eq = Join[
Flatten[Table[(λh uxx[xcol[[i]], ycol[[j]],
zcol[[k]]] + λc uyy[xcol[[i]], ycol[[j]],
zcol[[k]]] + λz uzz[xcol[[i]], ycol[[j]],
zcol[[k]]]) == 0, {i, M}, {j, M}, {k, M}]],
Flatten[Table[
u1[xcol[[i]], ycol[[j]], zcol[[k]]] -
u2[xcol[[i]], ycol[[j]], zcol[[k]]] == 0, {i, M}, {j, M}, {k,
M}]], Flatten[
Table[u1[xcol[[i]], ycol[[j]], zcol[[k]]] -
u3[xcol[[i]], ycol[[j]], zcol[[k]]] == 0, {i, M}, {j, M}, {k,
M}]], Flatten[
Table[ux[1, ycol[[i]], zcol[[j]]] == 0, {i, M}, {j, M}]],
Flatten[Table[uy[xcol[[i]], 1, zcol[[j]]] == 0, {i, M}, {j, M}]],
Flatten[Table[ux[0, ycol[[i]], zcol[[j]]] == 0, {i, M}, {j, M}]],
Flatten[Table[uy[xcol[[i]], 0, zcol[[j]]] == 0, {i, M}, {j, M}]],
Flatten[Table[
uz[xcol[[i]], ycol[[j]], 1] +
rh (u3[xcol[[i]], ycol[[j]], 1] - th[xcol[[i]], ycol[[j]]]) ==
0, {i, M}, {j, M}]],
Flatten[Table[
uz[xcol[[i]], ycol[[j]], 0] -
rc (u3[xcol[[i]], ycol[[j]], 0] - tc[xcol[[i]], ycol[[j]]]) ==
0, {i, M}, {j, M}]],
Flatten[Table[
thx[xcol[[i]], ycol[[j]]] -
bh (u3[xcol[[i]], ycol[[j]], 1] - th[xcol[[i]], ycol[[j]]]) ==
0, {i, M}, {j, M}]],
Flatten[Table[
tcy[xcol[[i]], ycol[[j]]] -
bc (u3[xcol[[i]], ycol[[j]], 0] - tc[xcol[[i]], ycol[[j]]]) ==
0, {i, M}, {j, M}]], Table[th[0, ycol[[i]]] == 1., {i, M}],
Table[tc[xcol[[i]], 0] == 0., {i, M}]];
var = Join[Flatten[U1], Flatten[U2], Flatten[U3], Flatten[B1],
Flatten[B2], Flatten[B3], Flatten[G1], Flatten[G2], Flatten[G3],
Flatten[G4], Flatten[G5], H1, H2];


Solution and visualization

{v, mat} = CoefficientArrays[eq, var];

sol1 = LinearSolve[mat, -v];

rul = Table[var[[i]] -> sol1[[i]], {i, Length[var]}];

{Plot3D[Evaluate[tc[x, y] /. rul], {x, 0, 1}, {y, 0, 1},
ColorFunction -> "Rainbow", MeshStyle -> White,
PlotLegends -> Automatic, PlotLabel -> \[Theta]c],
Plot3D[Evaluate[th[x, y] /. rul], {x, 0, 1}, {y, 0, 1},
ColorFunction -> "Rainbow", MeshStyle -> White,
PlotLegends -> Automatic, PlotLabel -> \[Theta]h],
Table[Plot3D[Evaluate[u1[x, y, z] /. rul], {x, 0, 1}, {y, 0, 1},
ColorFunction -> "Rainbow", MeshStyle -> White,
PlotLegends -> Automatic, PlotLabel -> \[Theta]w[z]], {z, 0,
1, .2}]}


We can compare this code with code developed by bbgodfrey above. But we need to change parameters as well

DSolveValue[{D[θc[y], y] + b (θc[y] - 1) ==
0, θc[0] == 0}, θc[y], y] // Simplify;
a00 = Simplify[Integrate[%, {y, 0, 1}]];
an0 = Simplify[Integrate[%% 2 Cos[n π y], {y, 0, 1}],
Assumptions -> n ∈ Integers];
DSolveValue[{D[θc[y], y] + b (θc[y] - Cos[m Pi y]) ==
0, θc[0] == 0}, θc[y], y] // Simplify;
a0m = Simplify[Integrate[%, {y, 0, 1}],
Assumptions -> m ∈ Integers];
amm = Simplify[Integrate[%% 2 Cos[m π y], {y, 0, 1}],
Assumptions -> m ∈ Integers];
anm = FullSimplify[Integrate[%%% 2 Cos[n π y], {y, 0, 1}],
Assumptions -> (m | n) ∈ Integers];
a[nn_?IntegerQ, mm_?IntegerQ] :=
Which[nn == 0 && mm == 0, a00, mm == 0, an0, nn == 0, a0m, nn == mm,
amm, True, anm] /. {n -> nn, m -> mm};

λh D[θw[x, y, z], x,
x] + λc D[θw[x, y, z], y,
y] + λz D[θw[x, y, z], z, z];
Simplify[(% /. θw ->
Function[{x, y, z},
Cos[nh Pi x] Cos[nc Pi y] θwz[z]])/(Cos[nh Pi x] Cos[
nc Pi y])] /. π^2 (nc^2 λc + nh^2 λh) ->
k[nh, nc]^2 λz;
Flatten@DSolveValue[% == 0, θwz[z], z] /. {C[1] -> c1[nh, nc],
C[2] -> c2[nh, nc]};

sz = c2[nh, nc] Sinh[k[nh, nc] z]/Cosh[k[nh, nc]] +
c1[nh, nc] Sinh[k[nh, nc] (1 - z)]/Sinh[k[nh, nc]];
sθc[nh_?IntegerQ, nc_?IntegerQ] :=
Sum[a[nc, m] c1[nh, m], {m, 0, maxc}] /. b -> bc;
(D[sz, z] == rc (sz - sθc[nh, nc])) /. z -> 0;
Solve[%, c2[nh, nc]] // Flatten // Apart;
sz1 = Simplify[sz /. %] // Apart;

szθh[nh_?IntegerQ, nc_?IntegerQ] := Evaluate[sz1 /. z -> 1];
sθh[nh_?IntegerQ, nc_?IntegerQ] :=
Evaluate[Sum[a[nh, m] szθh[m, nc], {m, 0, maxh}]];
eq = Simplify[(D[sz1, z] + rh (sz1 - sθh[nh, nc])) /. z -> 1] -
rh (DiscreteDelta[nh] - a[nh, 0]) DiscreteDelta[nc];

maxh = 3; maxc = 3; λh = 1; λc = 1; λz = 1; \
bh = 1; bc = 1; rh = 1; rc = 1;
ks = Flatten@
Table[k[nh, nc] ->
Sqrt[π^2 (nc^2 λc +
nh^2 λh)/λz], {nh, 0, maxh}, {nc, 0, maxc}];
eql = N[Collect[
Flatten@Table[
eq /. Sinh[k[0, 0]] -> k[0, 0], {nh, 0, maxh}, {nc, 0,
maxc}] /. b -> bh, _c1, Simplify] /. ks] /.
c1[z1_, z2_] :> Rationalize[c1[z1, z2]];
Union@Cases[eql, _c1, Infinity];
coef = NSolve[Thread[eql == 0], %] // Flatten;
sol = Total@
Simplify[
Flatten@Table[
sz1 Cos[nh Pi x] Cos[nc Pi y] /.
Sinh[z k[0, 0]] -> z k[0, 0], {nh, 0, maxh}, {nc, 0, maxc}],
Trig -> False] /. ks /. %;
(*{c1[0,0]->0.3788,c1[0,1]->-0.0234913,c1[0,2]->-0.00123552,c1[0,3]->-\
0.00109202,c1[1,0]->0.00168554,c1[1,1]->-0.0000775391,c1[1,2]->-5.\
40917*10^-6,c1[1,3]->-4.63996*10^-6,c1[2,0]->4.19045*10^-6,c1[2,1]->-\
1.24251*10^-7,c1[2,2]->-1.17696*10^-8,c1[2,3]->-1.02576*10^-8,c1[3,0]->\
1.65131*10^-7,c1[3,1]->-3.41814*10^-9,c1[3,2]->3.86348*10^-10,c1[3,3]->\
-3.48432*10^-10}*)


Visualization

pl1 = Table[
Plot3D[sol, {x, 0, 1}, {y, 0, 1}, PlotLegends -> Automatic,
AxesLabel -> {x, y}, Mesh -> None, ColorFunction -> "Rainbow",
PlotLabel -> Row[{"z = ", z}]], {z, 0, 1, .2}]


Numerical table for comparison on the grid

tab1 =
Table[{x, y, z, sol}, {x, 0, 1, .2}, {y, 0, 1, .2}, {z, 0, 1, .2}]

(*Out[]= {{{{0., 0., 0., 0.354583}, {0., 0., 0.2, 0.416427}, {0., 0.,
0.4, 0.472633}, {0., 0., 0.6, 0.527367}, {0., 0., 0.8,
0.583573}, {0., 0., 1., 0.645417}}, {{0., 0.2, 0., 0.361377}, {0.,
0.2, 0.2, 0.419296}, {0., 0.2, 0.4, 0.474032}, {0., 0.2, 0.6,
0.528107}, {0., 0.2, 0.8, 0.584005}, {0., 0.2, 1.,
0.645725}}, {{0., 0.4, 0., 0.375098}, {0., 0.4, 0.2,
0.426074}, {0., 0.4, 0.4, 0.47755}, {0., 0.4, 0.6, 0.530014}, {0.,
0.4, 0.8, 0.585128}, {0., 0.4, 1., 0.646529}}, {{0., 0.6, 0.,
0.387889}, {0., 0.6, 0.2, 0.433593}, {0., 0.6, 0.4,
0.481704}, {0., 0.6, 0.6, 0.532322}, {0., 0.6, 0.8,
0.586504}, {0., 0.6, 1., 0.647516}}, {{0., 0.8, 0.,
0.398835}, {0., 0.8, 0.2, 0.439582}, {0., 0.8, 0.4,
0.484998}, {0., 0.8, 0.6, 0.534165}, {0., 0.8, 0.8,
0.587608}, {0., 0.8, 1., 0.648311}}, {{0., 1., 0., 0.403914}, {0.,
1., 0.2, 0.441963}, {0., 1., 0.4, 0.486258}, {0., 1., 0.6,
0.534865}, {0., 1., 0.8, 0.588028}, {0., 1., 1.,
0.648613}}}, {{{0.2, 0., 0., 0.354275}, {0.2, 0., 0.2,
0.415995}, {0.2, 0., 0.4, 0.471893}, {0.2, 0., 0.6,
0.525968}, {0.2, 0., 0.8, 0.580704}, {0.2, 0., 1.,
0.638623}}, {{0.2, 0.2, 0., 0.361064}, {0.2, 0.2, 0.2,
0.418862}, {0.2, 0.2, 0.4, 0.473291}, {0.2, 0.2, 0.6,
0.526709}, {0.2, 0.2, 0.8, 0.581138}, {0.2, 0.2, 1.,
0.638936}}, {{0.2, 0.4, 0., 0.374776}, {0.2, 0.4, 0.2,
0.425637}, {0.2, 0.4, 0.4, 0.476809}, {0.2, 0.4, 0.6,
0.528616}, {0.2, 0.4, 0.8, 0.582265}, {0.2, 0.4, 1.,
0.639752}}, {{0.2, 0.6, 0., 0.38756}, {0.2, 0.6, 0.2,
0.433153}, {0.2, 0.6, 0.4, 0.480962}, {0.2, 0.6, 0.6,
0.530926}, {0.2, 0.6, 0.8, 0.583645}, {0.2, 0.6, 1.,
0.640754}}, {{0.2, 0.8, 0., 0.398498}, {0.2, 0.8, 0.2,
0.439139}, {0.2, 0.8, 0.4, 0.484255}, {0.2, 0.8, 0.6,
0.532769}, {0.2, 0.8, 0.8, 0.584753}, {0.2, 0.8, 1.,
0.641561}}, {{0.2, 1., 0., 0.403574}, {0.2, 1., 0.2,
0.441519}, {0.2, 1., 0.4, 0.485515}, {0.2, 1., 0.6,
0.53347}, {0.2, 1., 0.8, 0.585174}, {0.2, 1., 1.,
0.641868}}}, {{{0.4, 0., 0., 0.353471}, {0.4, 0., 0.2,
0.414872}, {0.4, 0., 0.4, 0.469986}, {0.4, 0., 0.6,
0.52245}, {0.4, 0., 0.8, 0.573926}, {0.4, 0., 1.,
0.624902}}, {{0.4, 0.2, 0., 0.360248}, {0.4, 0.2, 0.2,
0.417735}, {0.4, 0.2, 0.4, 0.471384}, {0.4, 0.2, 0.6,
0.523191}, {0.4, 0.2, 0.8, 0.574363}, {0.4, 0.2, 1.,
0.625224}}, {{0.4, 0.4, 0., 0.373936}, {0.4, 0.4, 0.2,
0.424502}, {0.4, 0.4, 0.4, 0.474899}, {0.4, 0.4, 0.6,
0.525101}, {0.4, 0.4, 0.8, 0.575498}, {0.4, 0.4, 1.,
0.626064}}, {{0.4, 0.6, 0., 0.3867}, {0.4, 0.6, 0.2,
0.432008}, {0.4, 0.6, 0.4, 0.47905}, {0.4, 0.6, 0.6,
0.527413}, {0.4, 0.6, 0.8, 0.576888}, {0.4, 0.6, 1.,
0.627096}}, {{0.4, 0.8, 0., 0.397621}, {0.4, 0.8, 0.2,
0.437988}, {0.4, 0.8, 0.4, 0.482342}, {0.4, 0.8, 0.6,
0.529259}, {0.4, 0.8, 0.8, 0.578005}, {0.4, 0.8, 1.,
0.627926}}, {{0.4, 1., 0., 0.402688}, {0.4, 1., 0.2,
0.440365}, {0.4, 1., 0.4, 0.483601}, {0.4, 1., 0.6,
0.52996}, {0.4, 1., 0.8, 0.578429}, {0.4, 1., 1.,
0.628242}}}, {{{0.6, 0., 0., 0.352484}, {0.6, 0., 0.2,
0.413496}, {0.6, 0., 0.4, 0.467678}, {0.6, 0., 0.6,
0.518296}, {0.6, 0., 0.8, 0.566407}, {0.6, 0., 1.,
0.612111}}, {{0.6, 0.2, 0., 0.359246}, {0.6, 0.2, 0.2,
0.416355}, {0.6, 0.2, 0.4, 0.469074}, {0.6, 0.2, 0.6,
0.519038}, {0.6, 0.2, 0.8, 0.566847}, {0.6, 0.2, 1.,
0.61244}}, {{0.6, 0.4, 0., 0.372904}, {0.6, 0.4, 0.2,
0.423112}, {0.6, 0.4, 0.4, 0.472587}, {0.6, 0.4, 0.6,
0.52095}, {0.6, 0.4, 0.8, 0.567992}, {0.6, 0.4, 1.,
0.6133}}, {{0.6, 0.6, 0., 0.385643}, {0.6, 0.6, 0.2,
0.430607}, {0.6, 0.6, 0.4, 0.476735}, {0.6, 0.6, 0.6,
0.523265}, {0.6, 0.6, 0.8, 0.569393}, {0.6, 0.6, 1.,
0.614357}}, {{0.6, 0.8, 0., 0.396543}, {0.6, 0.8, 0.2,
0.436578}, {0.6, 0.8, 0.4, 0.480024}, {0.6, 0.8, 0.6,
0.525113}, {0.6, 0.8, 0.8, 0.570518}, {0.6, 0.8, 1.,
0.615207}}, {{0.6, 1., 0., 0.4016}, {0.6, 1., 0.2,
0.438952}, {0.6, 1., 0.4, 0.481282}, {0.6, 1., 0.6,
0.525816}, {0.6, 1., 0.8, 0.570946}, {0.6, 1., 1.,
0.615531}}}, {{{0.8, 0., 0., 0.351689}, {0.8, 0., 0.2,
0.412392}, {0.8, 0., 0.4, 0.465835}, {0.8, 0., 0.6,
0.515002}, {0.8, 0., 0.8, 0.560418}, {0.8, 0., 1.,
0.601165}}, {{0.8, 0.2, 0., 0.358439}, {0.8, 0.2, 0.2,
0.415247}, {0.8, 0.2, 0.4, 0.467231}, {0.8, 0.2, 0.6,
0.515745}, {0.8, 0.2, 0.8, 0.560861}, {0.8, 0.2, 1.,
0.601502}}, {{0.8, 0.4, 0., 0.372074}, {0.8, 0.4, 0.2,
0.421995}, {0.8, 0.4, 0.4, 0.470741}, {0.8, 0.4, 0.6,
0.517658}, {0.8, 0.4, 0.8, 0.562012}, {0.8, 0.4, 1.,
0.602379}}, {{0.8, 0.6, 0., 0.384793}, {0.8, 0.6, 0.2,
0.429482}, {0.8, 0.6, 0.4, 0.474887}, {0.8, 0.6, 0.6,
0.519976}, {0.8, 0.6, 0.8, 0.563422}, {0.8, 0.6, 1.,
0.603457}}, {{0.8, 0.8, 0., 0.395675}, {0.8, 0.8, 0.2,
0.435447}, {0.8, 0.8, 0.4, 0.478174}, {0.8, 0.8, 0.6,
0.521826}, {0.8, 0.8, 0.8, 0.564553}, {0.8, 0.8, 1.,
0.604325}}, {{0.8, 1., 0., 0.400723}, {0.8, 1., 0.2,
0.437817}, {0.8, 1., 0.4, 0.479432}, {0.8, 1., 0.6,
0.522529}, {0.8, 1., 0.8, 0.564984}, {0.8, 1., 1.,
0.604656}}}, {{{1., 0., 0., 0.351387}, {1., 0., 0.2,
0.411972}, {1., 0., 0.4, 0.465135}, {1., 0., 0.6, 0.513742}, {1.,
0., 0.8, 0.558037}, {1., 0., 1., 0.596086}}, {{1., 0.2, 0.,
0.358132}, {1., 0.2, 0.2, 0.414826}, {1., 0.2, 0.4, 0.46653}, {1.,
0.2, 0.6, 0.514485}, {1., 0.2, 0.8, 0.558481}, {1., 0.2, 1.,
0.596426}}, {{1., 0.4, 0., 0.371758}, {1., 0.4, 0.2,
0.421571}, {1., 0.4, 0.4, 0.47004}, {1., 0.4, 0.6, 0.516399}, {1.,
0.4, 0.8, 0.559635}, {1., 0.4, 1., 0.597312}}, {{1., 0.6, 0.,
0.384469}, {1., 0.6, 0.2, 0.429054}, {1., 0.6, 0.4,
0.474184}, {1., 0.6, 0.6, 0.518718}, {1., 0.6, 0.8,
0.561048}, {1., 0.6, 1., 0.5984}}, {{1., 0.8, 0., 0.395344}, {1.,
0.8, 0.2, 0.435016}, {1., 0.8, 0.4, 0.477471}, {1., 0.8, 0.6,
0.520568}, {1., 0.8, 0.8, 0.562183}, {1., 0.8, 1.,
0.599277}}, {{1., 1., 0., 0.400389}, {1., 1., 0.2, 0.437385}, {1.,
1., 0.4, 0.478728}, {1., 1., 0.6, 0.521272}, {1., 1., 0.8,
0.562615}, {1., 1., 1., 0.599611}}}}*)


Results computed with our code for λh = 1; λc = 1; λz = 1; bh = 1; bc = 1; rh = 1; rc = 1;

Numerical table

tab2 =
Table[{x, y, z, u3[x, y, z] /. rul}, {x, 0, 1, .2}, {y, 0,
1, .2}, {z, 0, 1, .2}]

(*Out[]= {{{{0., 0., 0., 0.351913}, {0., 0., 0.2, 0.415671}, {0., 0.,
0.4, 0.472425}, {0., 0., 0.6, 0.527563}, {0., 0., 0.8,
0.584311}, {0., 0., 1., 0.648035}}, {{0., 0.2, 0., 0.362043}, {0.,
0.2, 0.2, 0.41954}, {0., 0.2, 0.4, 0.474254}, {0., 0.2, 0.6,
0.528512}, {0., 0.2, 0.8, 0.584869}, {0., 0.2, 1.,
0.648406}}, {{0., 0.4, 0., 0.374862}, {0., 0.4, 0.2,
0.426163}, {0., 0.4, 0.4, 0.477734}, {0., 0.4, 0.6,
0.530404}, {0., 0.4, 0.8, 0.585983}, {0., 0.4, 1.,
0.649196}}, {{0., 0.6, 0., 0.388098}, {0., 0.6, 0.2,
0.433759}, {0., 0.6, 0.4, 0.481916}, {0., 0.6, 0.6,
0.532732}, {0., 0.6, 0.8, 0.587369}, {0., 0.6, 1.,
0.650192}}, {{0., 0.8, 0., 0.398678}, {0., 0.8, 0.2,
0.439659}, {0., 0.8, 0.4, 0.485185}, {0., 0.8, 0.6,
0.534565}, {0., 0.8, 0.8, 0.588468}, {0., 0.8, 1.,
0.650978}}, {{0., 1., 0., 0.405092}, {0., 1., 0.2, 0.442489}, {0.,
1., 0.4, 0.486624}, {0., 1., 0.6, 0.535347}, {0., 1., 0.8,
0.588936}, {0., 1., 1., 0.651302}}}, {{{0.2, 0., 0.,
0.351575}, {0.2, 0., 0.2, 0.415109}, {0.2, 0., 0.4,
0.471481}, {0.2, 0., 0.6, 0.525731}, {0.2, 0., 0.8,
0.580452}, {0.2, 0., 1., 0.63791}}, {{0.2, 0.2, 0.,
0.361696}, {0.2, 0.2, 0.2, 0.418975}, {0.2, 0.2, 0.4,
0.47331}, {0.2, 0.2, 0.6, 0.526681}, {0.2, 0.2, 0.8,
0.581013}, {0.2, 0.2, 1., 0.638289}}, {{0.2, 0.4, 0.,
0.374505}, {0.2, 0.4, 0.2, 0.425594}, {0.2, 0.4, 0.4,
0.476789}, {0.2, 0.4, 0.6, 0.528573}, {0.2, 0.4, 0.8,
0.582132}, {0.2, 0.4, 1., 0.639099}}, {{0.2, 0.6, 0.,
0.387731}, {0.2, 0.6, 0.2, 0.433186}, {0.2, 0.6, 0.4,
0.48097}, {0.2, 0.6, 0.6, 0.530903}, {0.2, 0.6, 0.8,
0.583524}, {0.2, 0.6, 1., 0.640118}}, {{0.2, 0.8, 0.,
0.398303}, {0.2, 0.8, 0.2, 0.439083}, {0.2, 0.8, 0.4,
0.484238}, {0.2, 0.8, 0.6, 0.532737}, {0.2, 0.8, 0.8,
0.584628}, {0.2, 0.8, 1., 0.640923}}, {{0.2, 1., 0.,
0.404713}, {0.2, 1., 0.2, 0.441912}, {0.2, 1., 0.4,
0.485677}, {0.2, 1., 0.6, 0.53352}, {0.2, 1., 0.8,
0.585098}, {0.2, 1., 1., 0.641255}}}, {{{0.4, 0., 0.,
0.350792}, {0.4, 0., 0.2, 0.413998}, {0.4, 0., 0.4,
0.469592}, {0.4, 0., 0.6, 0.522253}, {0.4, 0., 0.8,
0.573834}, {0.4, 0., 1., 0.625097}}, {{0.4, 0.2, 0.,
0.360894}, {0.4, 0.2, 0.2, 0.417859}, {0.4, 0.2, 0.4,
0.47142}, {0.4, 0.2, 0.6, 0.523204}, {0.4, 0.2, 0.8,
0.574398}, {0.4, 0.2, 1., 0.625487}}, {{0.4, 0.4, 0.,
0.373682}, {0.4, 0.4, 0.2, 0.42447}, {0.4, 0.4, 0.4,
0.474897}, {0.4, 0.4, 0.6, 0.525099}, {0.4, 0.4, 0.8,
0.575526}, {0.4, 0.4, 1., 0.626319}}, {{0.4, 0.6, 0.,
0.386887}, {0.4, 0.6, 0.2, 0.432053}, {0.4, 0.6, 0.4,
0.479075}, {0.4, 0.6, 0.6, 0.527431}, {0.4, 0.6, 0.8,
0.576928}, {0.4, 0.6, 1., 0.627365}}, {{0.4, 0.8, 0.,
0.397442}, {0.4, 0.8, 0.2, 0.437944}, {0.4, 0.8, 0.4,
0.482342}, {0.4, 0.8, 0.6, 0.529268}, {0.4, 0.8, 0.8,
0.57804}, {0.4, 0.8, 1., 0.628191}}, {{0.4, 1., 0.,
0.403841}, {0.4, 1., 0.2, 0.440769}, {0.4, 1., 0.4,
0.48378}, {0.4, 1., 0.6, 0.530052}, {0.4, 1., 0.8,
0.578513}, {0.4, 1., 1., 0.628532}}}, {{{0.6, 0., 0.,
0.349794}, {0.6, 0., 0.2, 0.412617}, {0.6, 0., 0.4,
0.467266}, {0.6, 0., 0.6, 0.518075}, {0.6, 0., 0.8,
0.56624}, {0.6, 0., 1., 0.611867}}, {{0.6, 0.2, 0.,
0.359872}, {0.6, 0.2, 0.2, 0.416471}, {0.6, 0.2, 0.4,
0.469092}, {0.6, 0.2, 0.6, 0.519027}, {0.6, 0.2, 0.8,
0.566809}, {0.6, 0.2, 1., 0.612267}}, {{0.6, 0.4, 0.,
0.372633}, {0.6, 0.4, 0.2, 0.423072}, {0.6, 0.4, 0.4,
0.472566}, {0.6, 0.4, 0.6, 0.520924}, {0.6, 0.4, 0.8,
0.567945}, {0.6, 0.4, 1., 0.613119}}, {{0.6, 0.6, 0.,
0.385811}, {0.6, 0.6, 0.2, 0.430644}, {0.6, 0.6, 0.4,
0.476742}, {0.6, 0.6, 0.6, 0.523259}, {0.6, 0.6, 0.8,
0.569358}, {0.6, 0.6, 1., 0.614192}}, {{0.6, 0.8, 0.,
0.396346}, {0.6, 0.8, 0.2, 0.436526}, {0.6, 0.8, 0.4,
0.480006}, {0.6, 0.8, 0.6, 0.525098}, {0.6, 0.8, 0.8,
0.57048}, {0.6, 0.8, 1., 0.615039}}, {{0.6, 1., 0.,
0.40273}, {0.6, 1., 0.2, 0.439347}, {0.6, 1., 0.4,
0.481443}, {0.6, 1., 0.6, 0.525883}, {0.6, 1., 0.8,
0.570957}, {0.6, 1., 1., 0.615389}}}, {{{0.8, 0., 0.,
0.349011}, {0.8, 0., 0.2, 0.411519}, {0.8, 0., 0.4,
0.465435}, {0.8, 0., 0.6, 0.514807}, {0.8, 0., 0.8,
0.560343}, {0.8, 0., 1., 0.60129}}, {{0.8, 0.2, 0.,
0.359071}, {0.8, 0.2, 0.2, 0.415368}, {0.8, 0.2, 0.4,
0.46726}, {0.8, 0.2, 0.6, 0.51576}, {0.8, 0.2, 0.8,
0.560915}, {0.8, 0.2, 1., 0.601699}}, {{0.8, 0.4, 0.,
0.371809}, {0.8, 0.4, 0.2, 0.421961}, {0.8, 0.4, 0.4,
0.470731}, {0.8, 0.4, 0.6, 0.517659}, {0.8, 0.4, 0.8,
0.562058}, {0.8, 0.4, 1., 0.602568}}, {{0.8, 0.6, 0.,
0.384967}, {0.8, 0.6, 0.2, 0.429524}, {0.8, 0.6, 0.4,
0.474905}, {0.8, 0.6, 0.6, 0.519996}, {0.8, 0.6, 0.8,
0.563479}, {0.8, 0.6, 1., 0.603662}}, {{0.8, 0.8, 0.,
0.395485}, {0.8, 0.8, 0.2, 0.435399}, {0.8, 0.8, 0.4,
0.478168}, {0.8, 0.8, 0.6, 0.521837}, {0.8, 0.8, 0.8,
0.564608}, {0.8, 0.8, 1., 0.604526}}, {{0.8, 1., 0.,
0.401858}, {0.8, 1., 0.2, 0.438217}, {0.8, 1., 0.4,
0.479604}, {0.8, 1., 0.6, 0.522623}, {0.8, 1., 0.8,
0.565087}, {0.8, 1., 1., 0.604883}}}, {{{1., 0., 0.,
0.348701}, {1., 0., 0.2, 0.41105}, {1., 0., 0.4, 0.464655}, {1.,
0., 0.6, 0.513367}, {1., 0., 0.8, 0.557517}, {1., 0., 1.,
0.594879}}, {{1., 0.2, 0., 0.358753}, {1., 0.2, 0.2,
0.414897}, {1., 0.2, 0.4, 0.466479}, {1., 0.2, 0.6,
0.514321}, {1., 0.2, 0.8, 0.558091}, {1., 0.2, 1.,
0.595293}}, {{1., 0.4, 0., 0.371483}, {1., 0.4, 0.2,
0.421487}, {1., 0.4, 0.4, 0.46995}, {1., 0.4, 0.6, 0.51622}, {1.,
0.4, 0.8, 0.559237}, {1., 0.4, 1., 0.596173}}, {{1., 0.6, 0.,
0.384632}, {1., 0.6, 0.2, 0.429046}, {1., 0.6, 0.4,
0.474123}, {1., 0.6, 0.6, 0.518559}, {1., 0.6, 0.8,
0.560663}, {1., 0.6, 1., 0.597281}}, {{1., 0.8, 0.,
0.395142}, {1., 0.8, 0.2, 0.434919}, {1., 0.8, 0.4,
0.477385}, {1., 0.8, 0.6, 0.5204}, {1., 0.8, 0.8, 0.561795}, {1.,
0.8, 1., 0.598156}}, {{1., 1., 0., 0.401512}, {1., 1., 0.2,
0.437734}, {1., 1., 0.4, 0.478821}, {1., 1., 0.6, 0.521186}, {1.,
1., 0.8, 0.562276}, {1., 1., 1., 0.598518}}}}*)
{ListPlot[{Flatten[tab1, 2][[All, 4]], Flatten[tab2, 2][[All, 4]]},
PlotLegends -> {"tab1", "tab2"}, ImageSize -> Large,
PlotRange -> All],
ListPlot[Flatten[tab1, 2][[All, 4]] - Flatten[tab2, 2][[All, 4]],
PlotRange -> All]}


Note, that the difference of two methods is about $$2\times 10^{-3}$$.

• This works neatly. I did check with other configurations and materials to find solution in all cases. I see that you have used trigononetric basis functions with a combination of the wavelets. The pdf document uploaded by you on the wolfra community serves as a useful compendium. Jan 30, 2022 at 12:11
• @Avrana Yes, we can solve this problem with Euler wavelets as well. Jan 30, 2022 at 14:21

This is more of an extended comment than an answer, but it occurred to me that your solution is incomplete. You have a double $$Cos$$ series in $$m$$ and $$n$$, and unlike $$Sin$$ series you should need $$m=0$$ and $$n=0$$ terms.

You have computed your $$T_{mn}$$ series for $$(m, n)$$ going from $$1$$ to $$\infty$$ and it came out to be $$0$$. You need to add a $$T_{00}$$ term for $$(m, n)=0$$ and two more series.

Add a $$T_{m0}$$ series for $$n=0$$ and $$m$$ going from $$1$$ to $$\infty$$ and a $$T_{0n}$$ series for $$m=0$$ and n going from $$1$$ to $$\infty$$.

It takes all four pieces to make a complete solution. I have not tried this on your problem yet, so I don't know if all the pieces will come out to be zero or not, but this will give you something else to try. Your solution would not be correct without all four pieces anyway.

At the OP's request I will include my code, even though it doesn't work very well.

Clear["Global*"]
$Assumptions = n ∈ Integers && m ∈ Integers pde = D[T[x, y, z], x, x] + D[T[x, y, z], y, y] + D[T[x, y, z], z, z] == 0 T[x_, y_, z_] = X[x] Y[y] Z[z] pde = pde/T[x, y, z] // Expand x0eq = X''[x]/X[x] == 0 DSolve[x0eq, X[x], x] // Flatten X0 = X[x] /. % /. {C[1] -> c1, C[2] -> c2} xeq = X''[x]/X[x] == -α1^2 DSolve[xeq, X[x], x] // Flatten X1 = X[x] /. % /. {C[1] -> c3, C[2] -> c4} y0eq = Y''[y]/Y[y] == 0 DSolve[y0eq, Y[y], y] // Flatten Y0 = Y[y] /. % /. {C[1] -> c5, C[2] -> c6} yeq = Y''[y]/Y[y] == -β1^2 DSolve[yeq, Y[y], y] // Flatten Y1 = Y[y] /. % /. {C[1] -> c7, C[2] -> c8} z0eq = pde /. X''[x]/X[x] -> 0 /. Y''[y]/Y[y] -> 0 DSolve[z0eq, Z[z], z] // Flatten Z0 = Z[z] /. % /. {C[1] -> c9, C[2] -> c10} zeq = pde /. X''[x]/X[x] -> -α1^2 /. Y''[y]/Y[y] -> -β1^2 DSolve[zeq, Z[z], z] // Flatten Z1 = Z[z] /. % /. {C[1] -> c11, C[2] -> c12} // ExpToTrig // Collect[#, {Cosh[_], Sinh[_]}] & Z1 = % /. {c11 - c12 -> c11, c11 + c12 -> c12} T[x_, y_, z_] = X0 Y0 Z0 + X1 Y1 Z1 (D[T[x, y, z], x] /. x -> 0) == 0 c2 = 0; c4 = 0; T[x, y, z] c1 = 1 c3 = 1 (D[T[x, y, z], x] /. x -> L) == 0 α1 = (n π)/L (D[T[x, y, z], y] /. y -> 0) == 0 c6 = 0 c8 = 0 T[x, y, z] c5 = 1 c7 = 1 (D[T[x, y, z], y] /. y -> l) == 0 β1 = (m π)/l Tmn[x_, y_, z_] = T[x, y, z] /. {c9 -> 0, c10 -> 0} T00[x_, y_, z_] = T[x, y, z] /. n -> 0 /. m -> 0 T00[x_, y_, z_] = % /. c9 -> 0 /. c12 -> c1200 Tm0[x_, y_, z_] = T[x, y, z] /. n -> 0 Tm0[x_, y_, z_] = % /. {c10 -> 0, c9 -> 0, c11 -> c11m0, c12 -> c12m0} // PowerExpand T0n[x_, y_, z_] = T[x, y, z] /. m -> 0 // PowerExpand T0n[x_, y_, z_] = % /. {c9 -> 0, c10 -> 0, c11 -> c110n, c12 -> c120n} pdetcmn = D[tcmn[x, y], y] + (bc/l)*(tcmn[x, y] - Tmn[x, y, 0]) == 0 DSolve[pdetcmn, tcmn[x, y], {x, y}] // Flatten tcmn[x_, y_] = tcmn[x, y] /. % /. C[1][x] -> 0 pdetc00 = D[tc00[x, y], y] + (bc/l)*(tc00[x, y] - T00[x, y, 0]) == 0 DSolve[{pdetc00, tc00[x, 0] == tci}, tc00[x, y], {x, y}] // Flatten // Simplify tc00[x_, y_] = tc00[x, y] /. % pdetcm0 = D[tcm0[x, y], y] + (bc/l)*(tcm0[x, y] - Tm0[x, y, 0]) == 0 DSolve[pdetcm0, tcm0[x, y], {x, y}] // Flatten tcm0[x_, y_] = tcm0[x, y] /. % /. C[1][x] -> 0 pdetc0n = D[tc0n[x, y], y] + (bc/l)*(tc0n[x, y] - T0n[x, y, 0]) == 0 DSolve[pdetc0n, tc0n[x, y], {x, y}] // Flatten tc0n[x_, y_] = tc0n[x, y] /. % /. C[1][x] -> 0 pdethmn = D[thmn[x, y], x] + (bh/L)*(thmn[x, y] - Tmn[x, y, 0]) == 0 DSolve[pdethmn, thmn[x, y], {x, y}] // Flatten thmn[x_, y_] = thmn[x, y] /. % /. C[1][y] -> 0 pdeth00 = D[th00[x, y], x] + (bh/L)*(th00[x, y] - T00[x, y, 0]) == 0 DSolve[{pdeth00, th00[0, y] == thi}, th00[x, y], {x, y}] // Flatten th00[x_, y_] = th00[x, y] /. % pdethm0 = D[thm0[x, y], x] + (bh/L)*(thm0[x, y] - Tm0[x, y, 0]) == 0 DSolve[pdethm0, thm0[x, y], {x, y}] // Flatten thm0[x_, y_] = thm0[x, y] /. % /. C[1][y] -> 0 pdeth0n = D[th0n[x, y], x] + (bh/L)*(th0n[x, y] - T0n[x, y, 0]) == 0 DSolve[pdeth0n, th0n[x, y], {x, y}] // Flatten th0n[x_, y_] = th0n[x, y] /. % /. C[1][y] -> 0 bc100 = Simplify[(D[T00[x, y, z], z] /. z -> 0) == pc*(T00[x, y, 0] - tc00[x, y])] orth100 = Integrate[bc100[[1]], {y, 0, l}, {x, 0, L}] == Integrate[bc100[[2]], {y, 0, l}, {x, 0, L}] bc200 = Simplify[(D[T00[x, y, z], z] /. z -> w) == ph*(th00[x, y] - T00[x, y, w])] orth200 = Integrate[bc200[[1]], {y, 0, l}, {x, 0, L}] == Integrate[bc200[[2]], {y, 0, l}, {x, 0, L}] sol00 = Solve[{orth100, orth200}, {c10, c1200}] // Flatten // Simplify c10 = c10 /. sol00 c1200 = c1200 /. sol00 T00[x, y, z] tc00[x, y] th00[x, y] bc1m0 = Simplify[(D[Tm0[x, y, z], z] /. z -> 0) == pc*(Tm0[x, y, 0] - tcm0[x, y])] orth1m0 = Integrate[bc1m0[[1]]*Cos[(m*Pi*y)/l], {y, 0, l}, {x, 0, L}] == Integrate[bc1m0[[2]]*Cos[(m*Pi*y)/l], {y, 0, l}, {x, 0, L}] bc2m0 = Simplify[(D[Tm0[x, y, z], z] /. z -> w) == ph*(thm0[x, y] - Tm0[x, y, w])] orth2m0 = Integrate[bc2m0[[1]]*Cos[(m*Pi*y)/l], {y, 0, l}, {x, 0, L}] == Integrate[bc2m0[[2]]*Cos[(m*Pi*y)/l], {y, 0, l}, {x, 0, L}] solm0 = Solve[{orth1m0, orth2m0}, {c11m0, c12m0}] // Flatten // Simplify bc10n = (D[T0n[x, y, z], z] /. z -> 0) == pc*(T0n[x, y, 0] - tc0n[x, y]) orth10n = Integrate[bc10n[[1]]*Cos[(Pi*n*x)/L], {y, 0, l}, {x, 0, L}] == Integrate[bc10n[[2]]*Cos[(Pi*n*x)/L], {y, 0, l}, {x, 0, L}] bc20n = Simplify[(D[T0n[x, y, z], z] /. z -> w) == ph*(th0n[x, y] - T0n[x, y, w])] orth20n = Integrate[bc20n[[1]]*Cos[(Pi*n*x)/L], {y, 0, l}, {x, 0, L}] == Integrate[bc20n[[2]]*Cos[(Pi*n*x)/L], {y, 0, l}, {x, 0, L}] sol0n = Solve[{orth10n, orth20n}, {c110n, c120n}] // Flatten // Simplify bc1mn = (D[Tmn[x, y, z], z] /. z -> 0) == pc*(Tmn[x, y, 0] - tcmn[x, y]) orth1mn = Integrate[bc1mn[[1]]*Cos[(m*Pi*y)/l]*Cos[(Pi*n*x)/L], {y, 0, l}, {x, 0, L}] == Integrate[bc10n[[2]]*Cos[(m*Pi*y)/l]*Cos[(Pi*n*x)/L], {y, 0, l}, {x, 0, L}] bc2mn = Simplify[(D[Tmn[x, y, z], z] /. z -> w) == ph*(thmn[x, y] - Tmn[x, y, w])] orth2mn = Integrate[bc2mn[[1]]*Cos[(m*Pi*y)/l]*Cos[(Pi*n*x)/L], {y, 0, l}, {x, 0, L}] == Integrate[bc2mn[[2]]*Cos[(m*Pi*y)/l]*Cos[(Pi*n*x)/L], {y, 0, l}, {x, 0, L}] solmn = Solve[{orth1mn, orth2mn}, {c11, c12}] // Flatten // Simplify  All zeros except T00, and that solution does not satisfy the bc's. Have fun Update for new bc's This is too numerically unstable to get to work, but this is what I did. Clear["Global*"] pde = D[T[x, y, z], x, x] + D[T[x, y, z], y, y] + D[T[x, y, z], z, z] == 0$Assumptions = n ∈ Integers && m ∈ Integers && l > 0 && w > 0 && L > 0


Case 1

x = 0, T = thi

x = L, dT/dx = 0

y = 0, T = 0

y = l, dT/dy = 0 Use exponential in x, sinusoidal in y and z. Start with

T[x_, y_, z_] = (c1 + c2 x) (c10 z + c9) (c5 + c6 y) + (c3 Cosh[Sqrt[α1^2 + β1^2] x] +
c4 Sinh[Sqrt[α1^2 + β1^2] x]) (c7 Cos[α1 y] + c8 Sin[α1 y]) (c11 Sin[β1 z] + c12 Cos[β1 z])
T[0, y, z] == thi
(D[T[x, y, z], x] /. x -> L) == 0
c2 = 0
Solve[(c3 Sqrt[α1^2 + β1^2]Sinh[L Sqrt[α1^2 + β1^2]] +
c4 Sqrt[α1^2 + β1^2] Cosh[L Sqrt[α1^2 + β1^2]]) == 0, c4] // Flatten
c4 = c4 /. %
c3 = 1
c1 = 1


Manually expand the Tanh and incorporate the (constant) common denominator with the other constants

Simplify[Cosh[L*Sqrt[α1^2 + β1^2]]*Cosh[x*Sqrt[α1^2 + β1^2]] - Sinh[L*Sqrt[α1^2 + β1^2]]*Sinh[x*Sqrt[α1^2 + β1^2]]]
T[x_, y_, z_] = T[x, y, z] /. (Cosh[x Sqrt[α1^2 + β1^2]] -
Tanh[L Sqrt[α1^2 + β1^2]] Sinh[ x Sqrt[α1^2 + β1^2]]) -> %
T[x, 0, z] == 0
c5 = 0
c7 = 0
c6 = 1
c8 = 1

Simplify[D[T[x, y, z], y] /. y -> l] == 0
c10 = 0
c9 = 0
α1 = ((2 n + 1) π)/(2 l)


Set

β1 = ((2 m + 1) π)/(2 w)
T1[x_, y_, z_] = T[x, y, z]


Case 2

x = 0, T = 0

x = L, dT/dx = 0

y = 0, T = tci

y = l, dT/dy = 0

Use exponential in x, sinusoidal in y and z and flip the y and z terms

T2[x_, y_, z_] =
Sin[(π (2 n + 1) x)/(2 L)] (c112 Sin[(π (2 m + 1) z)/(2 w)] +
c122 Cos[(π (2 m + 1) z)/(2 w)]) Cosh[(l - y) Sqrt[(π^2 (2 n + 1)^2)/(4 L^2) + (π^2 (2 m + 1)^2)/(4 w^2)]]
T[x_, y_, z_] = T1[x, y, z] + T2[x, y, z]
pdeth = D[th[x, y], x] + (bh/L)*(th[x, y] - T[x, y, w]) == 0
DSolve[{pdeth, th[0, y] == thi}, th[x, y], {x, y}] //
Flatten // Simplify
th[x_, y_] = th[x, y] /. % // Simplify
pdetc = Simplify[D[tc[x, y], y] + (bc/l)*(tc[x, y] - T[x, y, 0]) == 0]
DSolve[{pdetc, tc[x, 0] == tci}, tc[x, y], {x, y}] //
Flatten // Simplify
tc[x_, y_] = tc[x, y] /. %
bc1 = T[0, y, z] == thi
bc2 = T[x, 0, z] == tci
bc3 = Simplify[(D[T[x, y, z], z] /. z -> 0) == pc*(T[x, y, 0] - tc[x, y])]
bc4 = Simplify[(D[T[x, y, z], z] /. z -> w) == ph*(th[x, y] - T[x, y, w])]
bc1eq = Simplify[Integrate[(bc1[[1]] - bc1[[2]])*Sin[(Pi*(2*n + 1)*y)/(2*l)]*Sin[(Pi*(2*m + 1)*z)/(2*w)], {z, 0, w}, {y, 0, l}] == 0]
bc2eq = Simplify[Integrate[(bc2[[1]] - bc2[[2]])*Sin[(Pi*(2*n + 1)*x)/(2*L)]*Sin[(Pi*(2*m + 1)*z)/(2*w)], {z, 0, w}, {x, 0, L}] == 0]
bc3eq = Integrate[bc3[[1]]*Sin[(Pi*(2*n + 1)*y)/(2*l)]*Sin[(Pi*(2*n + 1)*x)/(2*L)], {y, 0, l}, {x, 0, L}] == 0
bc4eq = Integrate[bc4[[1]]*Sin[(Pi*(2*n + 1)*y)/(2*l)]*Sin[(Pi*(2*n + 1)*x)/(2*L)], {y, 0, l}, {x, 0, L}] == 0
Solve[bc1eq, c12] // Flatten // Simplify
c12 = c12 /. %
Solve[bc2eq, c122] // Flatten // Simplify
c122 = c122 /. %
Solve[bc4eq, c112] // Flatten;
c112 = c112 /. %
Solve[bc3eq, c11] // Flatten;
c11 = c11 /. %
values = {L -> 1/40, l -> 1/40, w -> 3/1000, bh -> 433/1000,
bc -> 433/1000, ph -> 6524/100, pc -> 6524/100, thi -> 120, tci -> 30};
C11 = Table[c11 /. values, {m, 0, 10}, {n, 0, 10}] // N[#, 50] &
C11 = Re[C11]


To get rid of the small imaginary component. Chop wipes out the real part also.

C12 = Table[c12 /. values, {m, 0, 11}, {n, 0, 11}] // N[#, 50] &
C12 = Re[C12]
C112 = Table[c112 /. values, {m, 0, 11}, {n, 0, 11}] // N[#, 50] &
C112 = Re[C112]
C122 = Table[c122 /. values, {m, 0, 11}, {n, 0, 11}] // N[#, 50] &
C122 = Re[C122]


Put it all together

T[x_, y_, z_] := Sum[Sin[(Pi*(2*n + 1)*y)/(2*l)]*(C11[[m + 1,n + 1]]*Sin[(Pi*(2*m + 1)*z)/(2*w)] + C12[[m + 1,n + 1]]*Cos[(Pi*(2*m + 1)*z)/(2*w)])*
Cosh[(L - x)*Sqrt[(Pi^2*(2*n + 1)^2)/(4*l^2) + (Pi^2*(2*m + 1)^2)/(4*w^2)]] + Sin[(Pi*(2*n + 1)*x)/(2*L)]*
Cosh[(l - y)*Sqrt[(Pi^2*(2*n + 1)^2)/(4*L^2) + (Pi^2*(2*m + 1)^2)/(4*w^2)]]*(C112[[m + 1,n + 1]]*Sin[(Pi*(2*m + 1)*z)/(2*w)] +
C122[[m + 1,n + 1]]*Cos[(Pi*(2*m + 1)*z)/(2*w)]), {m, 0, 10}, {n, 0, 10}]


It took my computer days to compute all this and the values are way off. m,n of 10,10 are not enough terms, but I am not going any further. The values are still changing dramatically from m,n 9,10 to 10,10. Maybe the solution is wrong, or 50 decimals places is not enough, or it will take many more terms and many more days to even test the solution properly. Maybe your computer can do it faster, but my computer is 4 Ghz Intel i7 processor with 32 GB ram, so it is not a slow computer. Good luck.

• Yes, I did it separately. I am still examining the problem. I haven't gotten errors, but the bc's don't match. I will try some more. Jul 29, 2020 at 5:36
• My equations are so long that Mathematica runs all night and then runs out of memory in computing the c's. I will have to compute them with the numeric values applied before hand and then for each value of m, n separately. I am still working on it, though I still don't know if my answer will work yet. Aug 3, 2020 at 16:59
• Hopefully I will have something to post tomorrow, my time. Aug 6, 2020 at 7:30
• I am redoing it from a different perspective, but it will take a long time. Aug 9, 2020 at 5:52
• Happy new year to you too. Jan 1, 2021 at 8:27