# NDSolve "Coupled differential equation"

I am trying to solve the given coupled differential equation using NDsolve.

$$\begin{array}{c} n''(r)+\frac{n'(r) T'(r)}{T(r)}+\frac{n'(r)}{r}+\frac{n(r) T''(r)}{T(r)}-\frac{n(r) T'(r)^2}{T(r)^2}+\frac{n(r) T'(r)}{r T(r)}-n(r)+e^{-r^2}=0\\ T(r) n''(r)+2 n'(r) T'(r)+\frac{T(r) n'(r)}{r}+n(r) T''(r)+\frac{n(r) T'(r)}{r}-n(r)+e^{-r^2}=0 \end{array}$$

eq1 = Div[Grad[n[r] , {r, h, z}, "Cylindrical"], {r, h, z}, "Cylindrical"] + Div[(n[r]/T[r])*Grad[T[r], {r, h, z}, "Cylindrical"], {r, h, z}, "Cylindrical"] - n[r] + Exp[-(r^2)] == 0 ;

eq2 = Div[T[r]*Grad[n[r], {r, h, z}, "Cylindrical"], {r, h, z}, "Cylindrical"] + Div[n[r]*Grad[T[r] , {r, h, z}, "Cylindrical"], {r, h, z}, "Cylindrical"] - n[r] + Exp[-(r^2)] == 0 ;

initialConditions = {n'[50] == 0, n'[0] == 0, T'[50] == 0, T'[0] == 0};

solution = NDSolve[{eq1, eq2, initialConditions}, {n, T}, {r, 0, 50}];

Plot[Evaluate[{n[r], T[r]} /. %], {r, 0, 50}, PlotLegends -> {"n[r]", "T[r]"}]


Mathematica returns an error - NDSolve::bvdae: Differential-algebraic equations must be given as initial value problems

I have two questions regarding it,
Q1: Can we solve such equations accurately with NDSolve?
Q2: How do I get initial values for the two quantities, n and T?

• In V13.3.0, I get two messages, and the other one, which is printed first, seems more important. Did you not get that one? Commented Jul 26, 2023 at 14:36
• If you get the other message, begin to address it by evaluating Solve[{eq1, eq2}, {n''[r], T''[r]}]. It suggests there's an error in your ODEs. And if no error, then the ODEs are inconsistent and there's no solution at all. (The first error is the source of the second in my trial, which is why you should address the first error first.) Commented Jul 26, 2023 at 14:40
• @Shalu Rani: (1) Would behove you to search site for "differential algebraic equations" using search window above and read some of the threads, (2) Since cannot find explicit solutions for higher derivatives, may solve as IVP, (3) The DEs have $r$ and $T(r)$ in denominators so problematic at $r=0$, (4) Even if I set up as IVP with initial $r\neq 0$ again run into inconsistencies: Try solving for initial $T''[r_0]$ and $n''[r_0$] and cancelations of $T''[r_0]$ or $n''[r_0]$ occur likely causing the NDSolve::ivres message.
– josh
Commented Jul 26, 2023 at 15:46
• @Shalu Rani: Took the liberty to edit your post with the DEs in Latex.
– josh
Commented Jul 26, 2023 at 16:00
• @josh The DEs are inconsistent (the coefficients of $n''$ and $T''$ are proportional). You cannot fix that by adjusting ICs. The only hope is that the OP made a mistake in coding the DEs. Commented Jul 26, 2023 at 17:57

The pair of ODEs can be solved as follows. Begin by eliminating the second-derivative terms from eq2.

eq3 = Simplify[T[r] Subtract @@ eq1 - Subtract @@ eq2] == 0

(* E^-r^2 (-1 + T[r] - E^r^2 n'[r] T'[r]) + n[r] (1 - T[r] - T'[r]^2/T[r]) == 0 *)


Because n[r] and its derivatives enter the ODEs at most linearly, they can be eliminated to obtain an ODE for T[r] as follows.

Solve[{eq3, D[eq3, r], D[eq3, r, r]}, {n'''[r], n''[r], n'[r]}] //
Simplify // Flatten;
Simplify[{eq2, D[eq2, r]} /. %];
eqs = Equal @@@ (Simplify@Solve[%, {T'''[r], n[r]}] // Flatten)

(* {T’’’[r] == (-4 r^2 (-1 + r^2) T[r]^5 + r T[r]^4 (12 r (-1 + r^2)
+ (-7 + 8 r^2) T’[r] + r (5 - 4 r^2) T’’[r]) +
r T’[r]^2 ((-2 + 4 r^2) T’[r] + 4 r T’[r]^2 + r (-1 + 2 T’’[r])) +
T[r] ((2 + 4 r^2 - 4 r^4) T’[r]^2 - 8 r^3 T’[r]^3 + 2 r^2 T’[r]^4 -
r^2 (-1 + T’’[r]) T’’[r] + T’[r] (r - 2 r^3 - (r + 2 r^3) T’’[r])) +
T[r]^3 (2 (1 - 5 r^2 + 2 r^4) T’[r]^2 +
r T’[r] (15 - 18 r^2 + (-1 + 10 r^2) T’’[r]) +
r^2 (-12 (-1 + r^2) + (-9 + 8 r^2) T’’[r] - T’’[r]^2)) +
T[r]^2 ((2 r - 8 r^3) T’[r]^3 + T’[r]^2 (-4 + 7 r^2 - 2 r^2 T’’[r]) +
r T’[r] (-9 + 12 r^2 + (2 - 8 r^2) T’’[r]) +
r^2 (4 (-1 + r^2) + (3 - 4 r^2) T’’[r] +
2 T’’[r]^2)))/(r^2 (-1 + T[r])^2 T[r] (2 r T[r] - T’[r])),
n[r] == (E^-r^2 T[r] (r + r T[r]^2 + (1 - 2 r^2) T’[r] -
2 r T’[r]^2 + T[r] ((-1 + 2 r^2) T’[r] + r (-2 + T’’[r])) -
r T’’[r]))/(r T[r]^3 - r T’[r]^2 + T[r]^2 (-T’[r] + r (-2 + T’’[r])) +
T[r] (r + T’[r] - r T’[r]^2 - r T’’[r]))} *)


An important point is that it is a third-order system, so only three boundary conditions can be specified, not four as in the question. In the absence of suitable boundary conditions, I offer the following solution with made-up initial conditions:

NDSolveValue[{eqs, T[1] == 2, T'[1] == 0, T''[1] == 0}, {T[r], n[r]}, {r, 1, 3}];
Plot[%, {r, 1, 3}, AxesLabel -> {r, "T,n"}, LabelStyle -> {12, Bold, Black}]


Solution without eqs

An even simpler approach, at least for the sample solution, is

NDSolveValue[{eq2, eq3, T[1] == 2, T'[1] == 0, n'[1] == -.7},
{T[r], n[r]}, {r, 1, 3}];


which produces an answer very close to that plotted above.

• @bbjodfrey: Thanks for that interesting approach. I back-substituted your solutions into the original system over the range {r,1,3} and the error was on the order of $10^{-6}$, and no doubt improving with greater working precision.
– josh
Commented Jul 28, 2023 at 10:39
• @josh, after submitting my answer above, a ran the calculation above over a wider range and found that the simpler approach, "solution without eqs" was more robust. By the way, solving the boundary-value problem given in the question (of course, without one boundary condition) probably is possible, although challenging due to the strong nonlinearity of the ODEs. Commented Jul 30, 2023 at 1:04

This seems an overdetermined system.You get

     res1 = Solve[eq1 , T''[r] ][[1]]


$$\left\{T''(r)\to \frac{e^{-r^2} \left(-e^{r^2} r T(r)^2 n''(r)-e^{r^2} r T(r) n'(r) T'(r)-e^{r^2} T(r)^2 n'(r)-e^{r^2} n(r) T(r) T'(r)+e^{r^2} r n(r) T'(r)^2+e^{r^2} r n(r) T(r)^2-r T(r)^2\right)}{r n(r) T(r)}\right\}$$

But there is no solution for n''

    Simplify[eq2 /. res1 ] // FullSimplify


$$n'(r) T'(r)+n(r) \left(\frac{T'(r)^2}{T(r)}+T(r)-1\right)=e^{-r^2} (T(r)-1)$$

Using this equation as a replacement for n', there remains an undefined n in the equation for T''. Conclusion: There is an error in the algebraic ansatz.