I have the following code which solves for $(1)$ (i.e. solves for $C_1$ and $C_2$)
$$ T(x,y,z)=\sum_{n,m=0}^{\infty}(C_1 e^{\gamma z}+C_2 e^{-\gamma z})\sin\bigg(\frac{\alpha_n x}{L}+\beta_n\bigg)\sin\bigg(\frac{\delta_m y}{l}+\theta_m\bigg)+T_a \tag 1 $$
T[x_, y_, z_] = (C1* E^(γ z) + C2 E^(- γ z))*Sin[(α x/L) + β]*Sin[(δ y/l) + θ] + Ta
tc[x_, y_] = E^(-NTUC* y/l)*{tci + (NTUC/l)*Integrate[E^(NTUC*s/l)*T[x, s, 0], {s, 0, y}]};
tc[x_, y_] = tc[x, y][[1]];
bc1 = (D[T[x, y, z], z] /. z -> 0) == pc (T[x, y, 0] - tc[x, y]);
ortheq1 = Integrate[bc1[[1]]*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}] == Integrate[bc1[[2]]*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}];
ortheq1 = ortheq1 // Simplify;
th[x_, y_] = E^(-NTUH*x/L)*{thi + (NTUH/L)*Integrate[E^(NTUH*s/L)*T[s, y, w], {s, 0, x}]};
th[x_, y_] = th[x, y][[1]];
bc2 = (D[T[x, y, z], z] /. z -> w) == ph (th[x, y] - T[x, y, w]);
ortheq2 = Integrate[bc2[[1]]*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}] == Integrate[bc2[[2]]*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}];
ortheq2 = ortheq2 // Simplify;
soln = Solve[{ortheq1, ortheq2}, {Subscript[C, 1], Subscript[C, 2]}];
CC1 = C1 /. soln[[1, 1]];
CC2 = C2 /. soln[[1, 2]];
expression1 := CC1;
c1[α_, β_, δ_, θ_, γ_] := Evaluate[expression1];
expression2 := CC2;
c2[α_, β_, δ_, θ_, γ_] := Evaluate[expression2];
The following relations hold, $\beta_n=\tan^{-1}(1.66\times10^4 \alpha_n)$ and $\delta_m=\tan^{-1}(8.33\times10^3 \theta_m)$
The n=0 values is $\alpha_0=0.01095$ and m=0 value is $\delta_0=0.01549$.
Subsequently from n=1 and m=1 it is known that $\alpha_n=n\pi$ and $\delta_m=m\pi$.
I want to build a function such that this summation can be automatically performed for the desired values of $n$ and $m$.
$T_a$ is added only once in the final $T(x,y,z)$. The rest of the constants along with the other functions I wish to calculate are given below:
L = 0.9; l = 1.8; w = 0.0003; NTUH = 17.394; NNTUC = 22.151; ph = 8.6; pc = 13.93;
γ = Sqrt[(α/L)^2 + (δ/l)^2];
thi=460;tci=300;Ta=380;
tc1[x_, y_] = E^(-NTUC* y/l)*{tci + (NTUC*/l)*Integrate[E^(NTUC* s/l)*(TWnet /. {y -> s, z -> 0}), {s, 0, y}]};
th1[x_, y_] = E^(-NTUH* x/L)*{thi + (NTUH/L)*Integrate[E^(NTUH* s/L)*(TWnet /. {x -> s, z -> w}), {s, 0, x}]};
Plot[tc1[x, l], {x, 0, L}]
Plot[th1[L, y], {y, 0, l}]
THotAvg = Integrate[th1[x, y]/l, {y, 0, l}];
TColdAvg = Integrate[tc1[x, y]/L, {x, 0, L}];
THotAvg /. x -> L
TColdAvg /. y -> l
Plot[THotAvg, {x, 0, L}]
Plot[TColdAvg, {y, 0, l}]
The term TWnet in the above code section is the final $T(x,y,z)$ function I desire. So if someone can make the final distribution as a function then terms like TWnet /. {y -> s, z -> 0} would be something like TWnet[x,s,0]
I hope I was able to clearly explain the requirements here.
NOTE: The first code section takes some time to execute
CONTEXTUAL INFORMATION
I am trying to solve $\nabla^2 T(x,y,z)=0$ defined on $x\in[0,L], y\in[0,l]$ and $z\in[0,w]$ subjected to the following boundary conditions:
$$k(\frac{\partial T(0,y,z)}{\partial x})=h_a(T(0,y,z)-T_a) \tag A$$
$$-k(\frac{\partial T(L,y,z)}{\partial x})=h_a(T(L,y,z)-T_a) \tag B$$
$$k(\frac{\partial T(x,0,z)}{\partial y})=h_a(T(x,0,z)-T_a)\tag C$$
$$-k(\frac{\partial T(x,l,z)}{\partial y})=h_a(T(x,l,z)-T_a) \tag D$$
$$\frac{\partial T(x,y,0)}{\partial z} = p_c\bigg(T(x,y,0)-e^{-\beta_c y/l}\left[t_{ci} + \frac{\beta_c}{l}\int_0^y e^{\beta_c s/l}T(x,s,0)ds\right]\bigg) \tag E$$
$$\frac{\partial T(x,y,w)}{\partial z} = p_h\bigg(e^{-\beta_h x/L}\left[t_{hi} + \frac{\beta_h}{L}\int_0^x e^{\beta_h s/L}T(x,s,w)ds\right]-T(x,y,w)\bigg) \tag F$$
Now under the conditions $A,B,C,D$, the solution form of the three-dimensional Laplacian is given by $(1)$
$\gamma=\sqrt{(\alpha/L)^2 + (\delta/L)^2}$ (Have not mentioned this explicitly in the original question, so I wrote it here).
In the first section of the code I apply the $z$ boundary conditions and use orthogonality to determine the constants $C_1, C_2$. I must mention here that I have already proven the orthogonality of $\sin\bigg(\frac{\alpha_n x}{L}+\beta_n\bigg)$ under the boundary conditions $A-D$ The values of $\alpha$ and $\beta$ are to be calculated using the following transcendental equation:
$$2\cot{\alpha}=\frac{k\alpha}{h_a L}-\frac{h_aL}{k\alpha}\tag G$$ $$\beta=\tan^{-1}(\frac{k \alpha}{h_a L})\tag H$$
Similar set of equation exists for $\delta$ and $\theta$
I only want solution in the limit of very small $h_a \rightarrow 0$ for which except the first $\alpha$ value all other values are $n\pi$. I have derived an expression to calculate the first value which is:
$$\alpha=\frac 1{\sqrt a} \left( 1+\frac{1}{3 a}-\frac{8}{45 a^2}+\frac{53}{630 a^3}+O\left(\frac{1}{a^4}\right)\right)$$
where $a=k/(2h_a L)$. But in any case, I have posted the numerical values in the original question.
Once I get the $T(x,y,z)$ my objective is to calculate $t_h$ and $t_c$ which are given by:
$$t_h=e^{-\beta_h x/L}\bigg(t_{hi} + \frac{\beta_h}{L}\int_0^x e^{\beta_h s/L}T(x,s,w)ds\bigg) \tag I$$
$$t_c=e^{-\beta_c y/l}\bigg(t_{ci} + \frac{\beta_c}{l}\int_0^y e^{\beta_c s/l}T(x,s,0)ds\bigg) \tag J$$
Origins of the b.c.$E,F$
Actual bc(s): $$\frac{\partial T(x,y,0)}{\partial z}=p_c (T(x,y,0)-t_c) \tag K$$ $$\frac{\partial T(x,y,w)}{\partial z}=p_h (t_h-T(x,y,w))\tag L$$
where $t_h,t_c$ are defined in the equation:
$$\frac{\partial t_c}{\partial y}+\frac{\beta_c}{l}(t_c-T(x,y,0))=0 \tag M$$ $$\frac{\partial t_h}{\partial x}+\frac{\beta_h}{L}(t_h-T(x,y,0))=0 \tag N$$
It is known that $t_h(x=0)=t_{hi}$ and $t_c(y=0)=t_{ci}$. I had solved $M,N$ using the method of integrating factors and used the given conditions to reach $I,J$ which were then substituted into the original b.c.(s) $K,L$ to reach $E,F$.
My attempt I have written the following script to carry out the summation:
γ[α_, δ_] = Sqrt[(α/L)^2 + (δ/l)^2];
L = 0.9; l = 1.8; w = 0.0003; NTUH = 17.394; NTUC = 22.151; ph = 8.6; pc = 13.93;
α0 = 0.01095439637; δ0 = 0.0154917784; β0 = 1.56532; θ0 = 1.56305;
thi = 460; tci = 300; Ta = 380;
V0 = ((c1[α0, β0, δ0, θ0, γ[α0, δ0]] *E^(γ[α0, δ0] *z) + c2[α0, β0, δ0, θ0, γ[α0, δ0]]* E^(-γ[α0, δ0] *z))*Sin[δ0*y/l + θ0] + Sum[(c1[α0, β0, m*\[Pi], 1.5708,γ[α0, m*\[Pi]]] *E^(γ[α0, m*\[Pi]] *z) + c2[α0, β0, m*\[Pi], 1.5708, γ[α0, m*\[Pi]]]*E^(-γ[α0, m*\[Pi]]* z))*Sin[m*\[Pi]*y/l + 1.5708], {m, 1, 5}])*Sin[α0*x/L + β0];
Vn = Sum[((c1[n*\[Pi], 1.5708, δ0, θ0, γ[n*\[Pi], δ0]] *E^(γ[n*\[Pi], δ0] *z) + c2[n*\[Pi], 1.5708, δ0, θ0, γ[n*\[Pi], δ0]]* E^(-γ[n*\[Pi], δ0]* z))*Sin[δ0*y/l + θ0] + Sum[(c1[n*\[Pi], 1.5708, m*\[Pi], 1.5708, γ[n*\[Pi], m*\[Pi]]] *E^(γ[n*\[Pi], m*\[Pi]] *z) + c2[n*\[Pi], 1.5708, m*\[Pi], 1.5708, γ[n*\[Pi], m*\[Pi]]]* E^(-γ[n*\[Pi], m*\[Pi]]* z))*Sin[m*\[Pi]*y/l + 1.5708], {m, 1, 5}])*Sin[n*\[Pi]*x/L +1.5708], {n, 1, 5}];
Vnet = V0 + Vn + Ta;
tcf[x_, y_] = E^(-NTUC* y/l)*{tci + (NTUC/l)*Integrate[E^(NTUC* s/l)*(Vnet /. {y -> s, z -> 0}), {s, 0, y}]};
thf[x_, y_] = E^(-NTUH* x/L)*{thi + (NTUH/L)*Integrate[E^(NTUH* s/L)*(Vnet /. {x -> s, z -> w}), {s, 0, x}]};
tcfavg = Integrate[tcf[x, y], {x, 0, L}]/L;
thfavg = Integrate[thf[x, y], {y, 0, l}]/l;
tcfavg /. y -> l // Chop
thfavg /. x -> L // Chop
The tcfavg and thfavg plots i get are also weird
And the outlet temperatures are
tcfavg /. y -> l // Chop
401.984
thfavg /. x -> L // Chop
344.348
{x,y,z}but on the right hand side not ay. $\endgroup$\beta. Apologies for this inconvenience. But I have observed that when this code is copied from here to MMA, it is interpreted correctly $\endgroup$