# Problem with NDSolve: the function diverges

I'm trying to numerically simulate heat flux in a 3d power cable model that i build. Essentially the equation is:

I can get an analytic solution from Mathematica easily and when i plot T(x) i get no trouble (putting the constraint T'(x)=0 and T'(l)=0):

second picture shows T'(x)=0

Because i want to add some nonlinear thermal behavior i had to use NDSlove, and first i try to solve the same equation that i wrote before to see if it works:

sol = NDSolve[{-T[x] + \[Tau] +
P[x]*(Subscript[\[Rho], t]) + (Subscript[\[Rho],
t])/(Subscript[\[Rho], c])*T''[x] == 0, T'[0] == 0,  T'[l] == 0}, T[x], {x, 0, l}]


and the numerical solution give me this:

i can't undertand what's the deal here and why it doesn't converge!!

Here complete Code:

P[x_] := 3*r*(\[Kappa] + \[Lambda]*((x - l/2)/(l/2))^2)^2
T0 = P[0]*Subscript[\[Rho], t]
\[Kappa] = 1597.14
\[Lambda] = 1788 - \[Kappa]
r = 0.0000133
\[Tau] = 20
Subscript[\[Rho], t] = 0.54877
Subscript[\[Rho], c] = 0.332502
l = 100000
\[Alpha] = 0
sol = NDSolve[{-h[x] + \[Tau] +
P[x]*(Subscript[\[Rho], t]) + (Subscript[\[Rho],
t])/(Subscript[\[Rho], c])*h''[x] == 0, h'[0] == 0,
h'[l] == 0}, h[x], {x, 0, l}]

• Works for me with my chosen parameters: i.stack.imgur.com/BP3lR.png -- You'll need to give complete code to get much help. – Michael E2 Nov 23 '18 at 14:39
• complete code added – Mattia Nov 23 '18 at 14:47
• Parameters not defined. How can we verify the solution? – Alex Trounev Nov 23 '18 at 15:14
• In your picture analytical solution does not satisfy boundary conditions T'[0] == 0, T'[l] == 0 – Alex Trounev Nov 23 '18 at 15:33
• yes it's satisfied, slope of T(x) goes to zero at 0 and l but you can't see it from the picture because T''[x] is zero at 13 and graphic is from 0 to 100000, this is the reason why you can't see it. anyway i need help! – Mattia Nov 23 '18 at 15:49

1)Increase the order of the equation to solve the Dirichlet problem; 2) map the solution to the interval (0,1) by replacing x->l*x; 2) use the method of the false transient. Then the system of equations and the solution are

P[x_] := 3*r*(\[Kappa] + \[Lambda]*(2*x - 1)^2)^2

\[Kappa] = 1597.14;
\[Lambda] = 1788 - \[Kappa];
r = 0.0000133;
\[Tau] = 20;
rt = 0.54877;
rc = 0.332502;
l = 100000;
k = rt/rc/l^2;

eq = {-h[x, t] + D[P[x], x]*rt + k*D[h[x, t], x, x] -
D[h[x, t], t] == 0};
ic = h[x, 0] == 0;
bc = {DirichletCondition[h[x, t] == 0, x == 0],
DirichletCondition[h[x, t] == 0, x == 1]};

sol = NDSolveValue[{eq, ic, bc}, h, {x, 0, 1}, {t, 0, 20}]
Plot3D[sol[x, t], {x, 0, 1}, {t, 1, 20}, PlotRange -> All,
Mesh -> None, ColorFunction -> Hue,
AxesLabel -> {"x", "t", "\!$$\*SubscriptBox[\(h$$, $$x$$]\)"}]

Plot[Table[sol[x, t], {t, 5, 20, 5}], {x, 0, 1}, PlotRange -> All,
AxesLabel -> {"x", "\!$$\*SubscriptBox[\(h$$, $$x$$]\)"}]


The solution of the problem reaches a stationary state already at t>1. Integrating over x, we find the solution to the original problem up to a constant

hS[x_] := NIntegrate[sol[x1, 20], {x1, 0, x}]

{Plot[hS[x], {x, 0, .01}], Plot[hS[x], {x, 1 - .01, 1}],
Plot[hS[x], {x, 0, 1},AxesLabel -> {"x", "T"}]}


We choose the integration constant = 90, then the solution is similar to the analytical one. Note that with numerical integration on the right-hand boundary, when x-> 1, a numerical instability arises, although the function sol[x1,20] is rather smooth.