# Help for solving NDSolve

I am trying to solve below 3rd order DE using Runge Kutta methed but we have problem at t=0,

ClassicalRungeKuttaCoefficients[4, prec_] :=
With[{amat = {{1/2}, {0, 1/2}, {0, 0, 1}},
bvec = {1/6, 1/3, 1/3, 1/6}, cvec = {1/2, 1/2, 1}},
N[{amat, bvec, cvec}, prec]]

σ = .07;
ρ = 1000; g = 10; σT = .3*10^-3;
lσ = Sqrt[σ/(ρ*g)]; ΔT = .02;
Cr = (σT*ΔT)/σ; h0 = 1*10^-3;
Amp = (3*lσ^2*Cr)/(2*h0^2); T0 = 300;
T = ΔT*Exp[-(t/(200*10^-6))^2];
θ = (T - T0)/ΔT;
solution1 =
NDSolve[{y''[t] + 1/t*y'[t] - y[t] ==
Amp*Integrate[θ/(y[t] + 1), {t, -.008, .008}],
y[-0.008] == 0, y[0.008] == 0}, y[t], {t, -.008, .008},
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4,
"Coefficients" -> ClassicalRungeKuttaCoefficients},
StartingStepSize -> 1/10000]


• Not sure I understand the question. Your equation contains the term y'[t]/t, so I'm not really sure why it's surprising to find an error at t=0. But the solution1 you get out seems pretty useable as far as I can see. What exactly is the problem at t=0 (since you get a solution at the end)? What are you wanting to do with the solution that you can't do now? Jul 22, 2017 at 11:44
• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Jul 23, 2017 at 4:15
• It looks like there's a contradiction between the boundary conditions at the top of your post, $h(\pm\infty)=0$, and the solution graph where $h(r)$ tends towards 1. Should the BCs be $h'(\pm\infty)=0$? Jul 29, 2017 at 13:34

Because the ODE is singular at t == 0, it seems prudent to solve the ODE with DSolve rather than NDSolve. Unfortunately, using DSolve directly fails. So, instead solve

s = DSolveValue[{y''[t] + 1/t*y'[t] - y[t] == c,
y[-.008] == 0, y[.008] == 0}, y[t], {t, -.008, .008}] // Chop
(* c (-1 + 0.999984 BesselJ[0, I t]) *)


where c is the value of the integral, determined as follows.

ss = s/c;
int[c_] := Amp NIntegrate[θ/(ss c + 1), {t, -.008, .008}]
sc = c /. FindRoot[int[c] == c, {c, 1}, Evaluated -> False] // Chop
(* -0.211214 *)


Thus the solution is ss sc

Plot[sc ss, {t, -.008, .008}, ImageSize -> Large,
AxesLabel -> {t, y}, LabelStyle -> Directive[Bold, Black, 12]]


• @MichaelE2 I share your discomfort with the answer. However, the only way to eliminate the discontinuity is to force the coefficient of BesselY to be zero, in which case the given boundary conditions cannot be satisfied. Jul 22, 2017 at 19:53
• Thanks , typo in range corrected(-..008 to .008), But after that y(t) does not goes to Asymptotically zero.We need only real solution. We can change boundary condition 0 to some finite no. Jul 22, 2017 at 20:35
• @GopalVerma That change makes a big difference. Jul 22, 2017 at 20:48
• But,It should goes to Asymptotically zero. Jul 22, 2017 at 21:32
• I guess my previous comment is academic now. I like the sc = FixedPoint[int, -1.] trick. BTW, you can get rid of the last imaginary I with FunctionExpand[s]. Jul 22, 2017 at 23:56