I have asked this question before, but this is my new attempt and so instead of cluttering the previous one, I am making a new post. I am trying to analytically solve a PDE ($\nabla^2 T(x,y)=0$) coupled with an ODE. The PDE is subjected to the following boundary conditions:

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

$$\frac{\partial T(x,0)}{\partial y}=\gamma \tag 2$$

$$\frac{\partial T(x,l)}{\partial y}=\beta (T(x,l)-t) \tag 3$$

where $t$ is governed by the ODE:

$$\frac{\partial t}{\partial x}+\alpha(t-T(x,l))=0 \tag 4$$

subjected to $t(x=0)=0$. I am trying separation of variables. I manipulated $(4)$ to express $t$ as $t=\alpha e^{-\alpha x}\Bigg(\int_0^x e^{\alpha s }T(s,l)\mathrm{d}s\Bigg)$ and substituted in $(3)$ while applying the 3rd b.c.

My attempt is (I must acknowledge Bill Watts here as I have used methods I learned from his answer on MMA SE):

pde = D[T[x, y], x, x] + D[T[x, y], y, y] == 0

(*product form*)
T[x_, y_] = X[x] Y[y]
pde/T[x, y] // Expand
xeq = X''[x]/X[x] == -a^2
DSolve[xeq, X[x], x] // Flatten
X[x_] = X[x] /. % /. {C[1] -> c1, C[2] -> c2}
yeq = Y''[y]/Y[y] == a^2
DSolve[yeq, Y[y], y] // Flatten
Y[y_] = (Y[y] /. % /. {C[1] -> c3, C[2] -> c4})

(*addition form*)
T[x_, y_] = Xp[x] + Yp[y]
xpeq = Xp''[x] == b
DSolve[xpeq, Xp[x], x] // Flatten
Xp[x_] = Xp[x] /. % /. {C[1] -> c5, C[2] -> c6}
ypeq = Yp''[y] + b == 0
DSolve[ypeq, Yp[y], y] // Flatten
Yp[y_] = Yp[y] /. % /. {C[1] -> 0, C[2] -> c7}

T[x_, y_] = X[x] Y[y] + Xp[x] + Yp[y]
pde // FullSimplify

(*Applying the first and second b.c.*)
(D[T[x, y], x] /. x -> 0) == 0
c6 = 0
c2 = 0
c1 = 1

(D[T[x, y], x] /. x -> L) == 0
b = 0
a = (n π)/L
$Assumptions = n ∈ Integers

(*Applying the third b.c.*)
(D[T[x, y], y] /. y -> 0) == γ
c4 = c4 /. Solve[Coefficient[%[[1]], Cos[(π n x)/L]] == 0, c4][[1]]
c7 = c7 /. Solve[c7 == γ, c7][[1]]
T[x, y] // Collect[#, c3] &

(*now splitting T[x,y] into two parts*)
T[x, y] /. n -> 0
T0[x_, y_] = 2 c3 + c5 + y γ /. c5 -> 0
Tn[x_, y_] = T[x, y] - T0[x, y] // Simplify

(*applying the fourth b.c. to each part individually and using orthogonality*)

bcfn0 = (D[T0[x, y], y] /. y -> l) == β (T0[x, l] - α E^(-α x) Integrate[E^(α s) T0[s, l], {s, 0, x}])
Integrate[bcfn0[[1]], {x, 0, L}] == Integrate[bcfn0[[2]], {x, 0, L}]
Solve[%, c3]
c3 = c3 /. %[[1]]

bcfn = (D[Tn[x, y], y] /. y -> l) == β (Tn[x, l] - α E^(-α x) Integrate[E^(α s) Tn[s, l], {s, 0, x}])
Solve[Integrate[bcfn[[1]]*Cos[(n*Pi*x)/L], {x, 0, L}] == Integrate[bcfn[[2]]*Cos[(n*Pi*x)/L], {x, 0, L}], c5];
c5 = c5 /. %[[1]];//FullSimplify
T0[x_, y_] = T0[x, y] // Simplify
Tn[x_, y_] = Tn[x, y] // Simplify

Now we declare some constants and compile the functions

α = 62.9/2;
β = 1807/390;
γ = 3091.67/390;
L = 0.060;
l = 0.003;

T[x_, y_, mm_] := T0[x, y] + Sum[Tn[x, y], {n, 1, mm}]

Plot[{Evaluate[T[x, 0, 10]], Evaluate[T[x, l/2, 10]], Evaluate[T[x, l, 10]]}, {x, 0, L}]

The plot results are extremely ambiguous. The solution is not even converging (as I increase the number of terms , the T value keeps increasing). I cannot figure out what I have done wrong. Since the $T$ results are completely out, I have not calculated $t$. I cannot figure out what I have done wrong.

enter image description here

  • $\begingroup$ Are you sure that the t[0] == 0 condition is correct? $\endgroup$
    – bbgodfrey
    Aug 27, 2020 at 2:28
  • $\begingroup$ If it helps, then this condition can be taken as $t=t_i$ Then the $t$ differential equation evaluates to $$t=\alpha e^{-\alpha x}\Bigg(\int_0^x e^{\alpha s }T(s,l)\mathrm{d}s+\frac{t_i}{\alpha}\Bigg)$$ $t_i$ is a constant whose value can be taken as 300 $\endgroup$
    – Avrana
    Aug 27, 2020 at 3:11

1 Answer 1


I can fix the problem of your solution increasing with increasing n, but that will not give you a solution. Rather than copy your entire solution I will start where I think the problem starts.

You have

T0[x_, y_] = 2 c3 + c5 + y γ /. c5 -> 0

Change that to

T0[x_, y_] = 2 c3 + c5 + y γ /. c3 -> 0
(*c5 + γ y*)


Tn[x_, y_] = T[x, y] - T0[x, y] // FullSimplify
(*2 c3 Cos[(π n x)/L] Cosh[(π n y)/L]*)

In your case you were carrying an extra constant term c5 with Tn which was being added for every term in your sum which is why your solution increased with each term. In my case I carry c5 as the constant term, but only with T0. The changes below will require changing solving for c5 with bcf0 and solving for c3 with bcfn.

This next problem I fear is insurmountable with the calculation of bcfn0.

bcfn0 = (D[T0[x, y], y] /. y -> l) == β (T0[x, l] - α E^(-α x) Integrate[E^(α s) T0[s, l], {s, 0, x}]) // FullSimplify
(*γ E^(α x) == β (c5 + γ l)*)

Examining this result, it is obvious that there is no constant value c5 can take to satisfy this equation.

Furthermore, with the new Tn the orthogonality equation will result in c3 = 0. This means that T will have no x dependence, which when you think about it makes sense, if T is to satisfy Laplace eq and have x derivatives equal to zero at both ends in the x direction.

If T has no x dependence, then its derivatives can also have no x dependence, but with the y derivative of T depending on t which has x dependence, we have a problem.

  • $\begingroup$ Hi Bill, thanks. I had actually tried with c3->0, to reach a Tn=0 scenario. But I had no idea why in my present solution I had no convergence. But this became clear from your explanation in this answer. I think the problem lies with the convective type b.c. which seems physically correct, but is causing problems. $\endgroup$
    – Avrana
    Aug 28, 2020 at 4:16
  • $\begingroup$ Also, have you had any success with the cross-flow problem ? $\endgroup$
    – Avrana
    Aug 28, 2020 at 4:19
  • 2
    $\begingroup$ I had the exact same problem with the 3D problem. Tn = 0 and T0 not satisfying the bc. $\endgroup$
    – Bill Watts
    Aug 28, 2020 at 4:24
  • $\begingroup$ @IndrasisMitra I do not believe that this problem has a solution. You might consider posting it to Mathematics.StackExchange, because it looks more like a math than a Mathematica problem. $\endgroup$
    – bbgodfrey
    Aug 28, 2020 at 18:34
  • 1
    $\begingroup$ I will look at it when I get the chance to to do so in detail. I think this paper does a better job on details, than the other one, but we'll see. As far as the other paper, if you don't have the Stehfest algorithm, you can find it at library.wolfram.com/infocenter/MathSource/2691 if you haven't discovered it already. $\endgroup$
    – Bill Watts
    Sep 1, 2020 at 0:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.