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[0] == A, y[Infinity] == 0}, y[x], x] 

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

  • 2
    $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Oct 27, 2021 at 9:15
  • $\begingroup$ 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[1] + x^((2 - D)/2) BesselY[1/2 (-2 + D), -I Sqrt[k] x] C[2]}}. After that how can I apply the boundary condition at infinity? $\endgroup$
    – New User
    Oct 27, 2021 at 9:27
  • $\begingroup$ Try converting the Bessel functions into modified Bessel functions and then recall that only BesselK vanishes at infinity. $\endgroup$
    – bbgodfrey
    Oct 27, 2021 at 9:31
  • $\begingroup$ 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. $\endgroup$
    – Nasser
    Oct 27, 2021 at 9:34
  • 2
    $\begingroup$ 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. $\endgroup$
    – Michael E2
    Oct 27, 2021 at 10:31

3 Answers 3


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

This is a guess

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]]

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)$$

  • $\begingroup$ It's brilliant! I believe this is what I am looking for. Thanks a lot! $\endgroup$
    – New User
    Oct 28, 2021 at 10:31

Can you please show your method?

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[0] == a, y[b] == 0};
sol = y[x] /. First@DSolve[{ode, ic}, y[x], x]

Mathematica graphics


 Assuming[d > 2, Simplify[sol]]
 (* 0 *)
  • $\begingroup$ 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? $\endgroup$
    – New User
    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.


ode = y''[x] + ((d - 1) Derivative[1][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[0] == 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

Limit[y[x], x -> 0]

Limit[y[x], x -> \[Infinity]]

Example plot

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

enter image description here


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.