# Analytical solution of second order linear differential equation with boundary at infinity

I am a new user to Mathematica and I would like to solve a simple second-order differential equation as follow:

$$y''[x]+\frac{(D-1)}{x}\times y'[x]=k\times y[x]$$,

where $$D$$ and $$k$$ are just two parameters and the boundary conditions are $$y[x=0]=A$$ and $$y[x=\infty]=0$$.

How can I get an analytical solution for the equation with the boundary condition at infinity?

My attempt for the question looks like:

sol = DSolve[{y''[x] + (D - 1)*y'[x]/x == k*y[x], y == A, y[Infinity] == 0}, y[x], x]


but the solver take the second input boundary condition as a 'True' argument.

• y[Infinity] == 0 is not a valid boundary condition. Try using DSolve[{y''[x] + (D - 1)*y'[x]/x == k*y[x]}, y[x], x]  and then applying the boundary conditions after the fact. Oct 27, 2021 at 9:15
• I used DSolve[{y''[x] + (D - 1)*y'[x]/x == k*y[x]}, y[x], x] and I get the following result : {{y[x] -> x^((2 - D)/2) BesselJ[1/2 (-2 + D), -I Sqrt[k] x] C + x^((2 - D)/2) BesselY[1/2 (-2 + D), -I Sqrt[k] x] C}}. After that how can I apply the boundary condition at infinity? Oct 27, 2021 at 9:27
• Try converting the Bessel functions into modified Bessel functions and then recall that only BesselK vanishes at infinity. Oct 27, 2021 at 9:31
• what is the value of D ? Btw, you should not use D but use d. I found the solution to be zero when d>2. otherwise, not defined. Oct 27, 2021 at 9:34
• Part of Nasser's point is that D is the partial derivative operator and a protected symbol. It's best practice to avoid single-letter capitals for your own variables. Oct 27, 2021 at 10:31

It is hard to get a fully automatic solution, but here is a human-assisted way:

This is a guess

Clear[ys]
ys[x_] := x^n  BesselK[n, Sqrt[k] x]


Does it work?

ys''[x] + (d - 1)*ys'[x]/x == k*ys[x] // FullSimplify


Seems so

(*Sqrt[k] (-2 + d + 2 n) x^n BesselK[-1 + n, Sqrt[k] x] == 0*)


But let us check the boundary conditions, taking into account $$n=\frac{2-d}{2}$$.

a=Assuming[k > 0 && n > 0,
Limit[ys[x], x -> 0]] /. {n -> (2 - d)/2} // Simplify


$$a=-2^{-\frac{d}{2}-1} d k^{\frac{d-2}{4}} \Gamma \left(-\tfrac{d}{2}\right)$$

Check the boundary condition at infinity

Assuming[k > 0 && n > 0, Limit[ys[x], x -> Infinity]]
(*0*)


Thus two boundary conditions are fulfilled and the solution reads

Full solution then $$y=\tfrac{A}{a} x^n K_n\left(\sqrt{k} x\right)$$

• It's brilliant! I believe this is what I am looking for. Thanks a lot! Oct 28, 2021 at 10:31

Too large to post as comment. will remove this if not useful. Basically, just solved the ode with b in place of infinity (with the idea of later taking the limit as b->infinity).

Then simplified the result with assumption d>2 which gives zero. No need to take limit. If d is not larger than 2, solution as given by Mathematica is not defined, since the solution to the ode has terms that look like

     0^(1/2 (-2 + d))


Which is not defined unless the power is positive (i.e. d>2) . And then it is zero. Here is the code

ClearAll[y, x, d, a, k, b];
ode = y''[x] + (d - 1)*y'[x]/x == k*y[x];
ic = {y == a, y[b] == 0};
sol = y[x] /. First@DSolve[{ode, ic}, y[x], x] Now

 Assuming[d > 2, Simplify[sol]]
(* 0 *)


• Hi, thanks for showing me the method. I think that the 0^(1/2 (-2 + d)) term might be due to the (d-1)/x term in the original equation, which makes the boundary condition of y[x=0]=a invalid. I can get a numerical solution by changing the boundary condition to y[x=0.0001]=a. So I would like to ask if it is also possible to use limit to indicate this first boundary condition (the second boundary condition with infinity x position remains the same) so that the analytical solution can be defined with 0<d<1? Oct 28, 2021 at 5:56

I don't know if there is a way to solve the ode for general d between 0 and 1, but it is possible to get a solution for a particular d. For example I will use d = 1/2.

Clear["Global*"]

ode = y''[x] + ((d - 1) Derivative[y][x])/x - k y[x] == 0 /. d -> 1/2


Using the finite bc's as a first pass is a good idea.

bc = {y == a, y[b] == 0}

DSolve[{ode, bc}, y[x], x] // Flatten // Simplify

y[x_] = y[x] /. %


We get some complexes that we don't want so go through some simplification routines.

y[x_] = Simplify[FunctionExpand[y[x]], k > 0]

y[x_] = PowerExpand[% // Expand]

\$Assumptions = k > 0


The above assumption is necessary to get the limit as b -> \[Infinity]

y[x_] = Limit[y[x], b -> \[Infinity]]
(*(2^(1/4) a k^(3/8) x^(3/4) BesselK[3/4, Sqrt[k] x])/Gamma[3/4]*)


Again, this solution is only valid for d =1/2.

Check the solution.

ode // FullSimplify
(*True*)

Limit[y[x], x -> 0]
(*a*)

Limit[y[x], x -> \[Infinity]]
(*0*)


Example plot

Plot[y[x] /. {a -> 1, k -> 1}, {x, 0, 10}, PlotRange -> All]
` 