How to solve these two differential equations

I have two differential equations and I try to use function DSolve to solve them together.

DSolve[{r^2*A''[r] + 2*r*A'[r] - 2 A[r] + lambda^2*r^2*A[r] +
2 lambda*r^2*dd'[r] == 0,
r*dd''[r] + 2 dd'[r] + 3 lambda^2*r*dd[r] + 2 lambda*r*A'[r] +
4 lambda*A[r] == 0}, {A, dd}, r]


There are two functions as A and dd. But it failed.

I try to use Maple software to solve it and it successes and needs little time.

The Maple code:

sol3 := r^2*diff(A(r), r, r) + 2*r*diff(A(r), r) - 2*A(r) + lambda^2*r^2*A(r) +2*lambda*r^2*diff(dd(r), r) = 0

sol4 := r*diff(dd(r), r, r) + 2*diff(dd(r), r) + 3*lambda^2*r*dd(r) + 2*lambda*r*diff(A(r), r) + 4*lambda*A(r) = 0

dsolve({sol3, sol4})


Therefore, the equations has solutions. I want to know whether Mathematica can solve it and how to do it. Thank you very much.

Edited to derive simpler expression for dd

The symbolic solution requested in the Question is obtained as follows.

eq1 = Expand[(r^2*A''[r] + 2*r*A'[r] - 2 A[r] + lambda^2*r^2*A[r] +
2 lambda*r^2*dd'[r])/r^2];
(* lambda^2 A[r] - 2 A[r]/r^2 + 2 A'[r]/r + 2 lambda dd'[r] + A''[r] *)
eq2 = Expand[(r*dd''[r] + 2 dd'[r] + 3 lambda^2*r*dd[r] +
2 lambda*r*A'[r] + 4 lambda*A[r])/r];
(* (4 lambda A[r])/r + 3 lambda^2 dd[r] + 2 lambda A'[r] + 2 dd'[r]/r + dd''[r] *)


Combine the two expressions to obtain an expression for dd in terms of A

eqdd = -Subtract @@ First[Eliminate[{eq1 == 0, eq2 == 0, D[eq1, r] == 0},
{dd'[r], dd''[r]}]]
(* -6*lambda^2*A[r]/r - 6*lambda^3*dd[r] - 3*lambda^2*A'[r] + 4*A''[r]/r + A'''[r] *)


and an expression for A alone and solve it.

eqA = Subtract @@ First[Eliminate[{D[eqdd, r] == 0, eq1 == 0}, dd'[r]]]
(* 3*lambda^4*r*A[r] - 4*A''[r]/r + 4*A'''[r] + r*A''''[r] *)

sA = DSolveValue[eqA == 0, A[r], r]
(*   (E^((-1)^(3/4) 3^(1/4) lambda r) (I + (-3)^(1/4) lambda r) C[1])/r^2
+ (E^((-3)^(1/4) lambda r) (-I + (-1)^(3/4) 3^(1/4) lambda r) C[2])/r^2
+ (E^(-(-1)^(3/4) 3^(1/4) lambda r) (-I + (-3)^(1/4) lambda r) C[3])/r^2
+ (E^(-(-3)^(1/4) lambda r) (I + (-1)^(3/4) 3^(1/4) lambda r) C[4])/r^2 *)


Finally, substitute A into eqdd and solve it.

sdd = Solve[Simplify[eqdd /. A -> Function[{r}, sA]] == 0, dd[r]][[1, 1]] // Values
(* -(1/(lambda r^3)) (E^((-1)^(3/4) 3^(1/4) lambda r) (I + (-3)^(1/4) lambda r) C[1] +
E^((-3)^(1/4) lambda r) (-I + (-1)^(3/4) 3^(1/4) lambda r) C[2] +
E^(-(-1)^(3/4) 3^(1/4) lambda r) (-I + (-3)^(1/4) lambda r) C[3] +
E^(-(-3)^(1/4) lambda r) (I + (-1)^(3/4) 3^(1/4) lambda r) C[4]) *)


as requested.

Addendum: Eliminating complex numbers from the solution.

It often is desirable not to have complex numbers in the solutions of real ODEs. This can be accomplished by converting exponentials to trigonometric and hyperbolic functions, performing a complex expansion, and replacing the existing constants of integration by constants that absorb the complex numbers.

sA1 = Collect[ComplexExpand@ExpToTrig@sA, {r, Cos[_] Cosh[_], Cos[_] Sinh[_],
Sin[_] Cosh[_], Sin[_] Sinh[_]}, Simplify];
Cases[sA1, Cos[_] Cosh[_] z_ | Sin[_] Cosh[_] z_ | Cos[_] Sinh[_] z_ |
Sin[_] Sinh[_] z_ -> z, Infinity, 4];
cf = Solve[Thread[{c1, c2, c3, c4} == %], {C[1], C[2], C[3], C[4]}] // Flatten;
sA2 = Collect[sA1 /. cf, {r, Cos[_] Cosh[_], Cos[_] Sinh[_], Sin[_] Cosh[_],
Sin[_] Sinh[_]}, Simplify]
(* (1/(r^2))(c1 Cos[(3^(1/4) lambda r)/Sqrt[2]] Cosh[(3^(1/4) lambda r)/ Sqrt[2]] +
c2 Cosh[(3^(1/4) lambda r)/Sqrt[2]] Sin[(3^(1/4) lambda r)/Sqrt[2]] +
c3 Cos[(3^(1/4) lambda r)/Sqrt[2]] Sinh[(3^(1/4) lambda r)/ Sqrt[2]] +
c4 Sin[(3^(1/4) lambda r)/Sqrt[2]] Sinh[(3^(1/4) lambda r)/Sqrt[2]]) +
(1/r)(-((3^(1/4) (c2 + c3) lambda Cos[(3^(1/4) lambda r)/Sqrt[2]] Cosh[(
3^(1/4) lambda r)/Sqrt[2]])/Sqrt[2]) +
(3^(1/4) (c1 - c4) lambda Cosh[(3^(1/4) lambda r)/Sqrt[2]] Sin[(
3^(1/4) lambda r)/Sqrt[2]])/Sqrt[2] -
(3^(1/4) (c1 + c4) lambda Cos[(3^(1/4) lambda r)/Sqrt[2]] Sinh[(
3^(1/4) lambda r)/Sqrt[2]])/Sqrt[2] -
(3^(1/4) (c2 - c3) lambda Sin[(3^(1/4) lambda r)/Sqrt[2]] Sinh[(
3^(1/4) lambda r)/Sqrt[2]])/Sqrt[2]) *)

sdd2 = Collect[ComplexExpand@ExpToTrig@sdd /. cf, {r, Cos[_] Cosh[_],
Cos[_] Sinh[_], Sin[_] Cosh[_], Sin[_] Sinh[_]}, Simplify]
(* (1/(r^2))((3^(1/4) (c2 + c3) Cos[(3^(1/4) lambda r)/Sqrt[2]] Cosh[(
3^(1/4) lambda r)/Sqrt[2]])/Sqrt[2] +
(3^(1/4) (-c1 + c4) Cosh[(3^(1/4) lambda r)/Sqrt[2]] Sin[(
3^(1/4) lambda r)/Sqrt[2]])/Sqrt[2] +
(3^(1/4) (c1 + c4) Cos[(3^(1/4) lambda r)/Sqrt[2]] Sinh[(
3^(1/4) lambda r)/Sqrt[2]])/Sqrt[2] +
(3^(1/4) (c2 - c3) Sin[(3^(1/4) lambda r)/Sqrt[2]] Sinh[(
3^(1/4) lambda r)/Sqrt[2]])/Sqrt[2]) +
(1/(r^3))(-((c1 Cos[(3^(1/4) lambda r)/Sqrt[2]] Cosh[(3^(1/4) lambda
r)/Sqrt[2]])/lambda) -
(c2 Cosh[(3^(1/4) lambda r)/Sqrt[2]] Sin[(3^(1/4) lambda
r)/Sqrt[2]])/lambda -
(c3 Cos[(3^(1/4) lambda r)/Sqrt[2]] Sinh[(3^(1/4) lambda
r)/Sqrt[2]])/lambda -
(c4 Sin[(3^(1/4) lambda r)/Sqrt[2]] Sinh[(3^(1/4) lambda
r)/Sqrt[2]])/lambda) *)


For convenience, the solutions can be plotted for one non-zero c at a time.

Plot[Evaluate[Replace[sA2, {{c1 -> 1}, {c2 -> 1}, {c3 -> 1}, {c4 -> 1}}, Infinity]
/. {c1 -> 0, c2 -> 0, c3 -> 0, c4 -> 0, lambda -> 1}], {r, 0, 5},
ImageSize -> Large, AxesLabel -> {r, A}, LabelStyle -> {15, Bold, Black},
PlotLegends -> Placed[{c1, c2, c3, c4}, {.3, .8}]]


Plot[Evaluate[Replace[sdd2, {{c1 -> 1}, {c2 -> 1}, {c3 -> 1}, {c4 -> 1}}, Infinity]
/. {c1 -> 0, c2 -> 0, c3 -> 0, c4 -> 0, lambda -> 1}], {r, 0, 5},
ImageSize -> Large, AxesLabel -> {r, dd}, LabelStyle -> {15, Bold, Black},
PlotLegends -> Placed[{c1, c2, c3, c4}, {.6, .8}]]


• +1. In the hands of a great chef . Commented Sep 14, 2020 at 5:35
• The question did not pose that A has to be differentiable four times. There are no conditions given to solve for all the constants. Problem is to solve this with minimal additional openness not to derive a polynomial over the Complexes. Great is, this uses only DSolve. The problem is totally underdetermined. Commented Sep 15, 2020 at 19:31
ClearAll[A, r, dd, sol3, sol4, lambda]
sol3 = r^2*D[A[r], r, r] + 2*r*D[A[r], r] - 2*A[r] +
lambda^2*r^2*A[r] + 2*lambda*r^2*D[dd[r], r]

sol4 = r*D[dd[r], r, r] + 2*D[dd[r], r] + 3*lambda^2*r*dd[r] +
2*lambda*r*D[A[r], r] + 4*lambda*A[r]
lambda = 0;

s = NDSolve[{sol3 == 0, sol4 == 0, A[1] == 1, A'[1] == 1, dd[1] == 1,
dd'[1] == 1}, {A, dd}, {r, 1, 2}]


Plot[{A[r], dd[r]} /. s, {r, 1, 2}, PlotRange -> All]


For

lambda = 0.5;


If the integration interval is extented towards zero the singularity at r==0 get important and influences, dominates the solution:

lambda = 5.0; s2 = NDSolve[{sol3 == 0, sol4 == 0, A1 == 1, A'1 == 1, dd1 == 1, dd'1 == 1}, {A, dd}, {r, MachinePrecision, 5}] Plot[{A[r], dd[r]} /. s2, {r, 0.2, 2}, PlotRange -> All, PlotStyle -> Automatic]

Mathematica has no DSolve solution since the details are not given in the question. The boundary conditions for this solution are taken by choice. They can be selected from a wide range. Same is true for lambda.

lambda = 5.0;
s2 = DSolve[{sol3 == 0, sol4 == 0, A[1] == 0, A'[1] == .1, dd[1] == 2,
dd'[1] == 1}, {A, dd}, {r, 0.1, 5}]
(*{{A -> Function[{r}, -(0.0333333/r^2) + 0.0333333 r],
dd -> Function[{r}, 3 - 1/r]}}*)

Plot[{A[r], dd[r]} /. s2, {r, 0.1, 5}, PlotRange -> All,
PlotStyle -> Automatic]


• Your DSolve answer is valid only for lambda = 0, as can be verified by Simplify[{sol3, sol4} /. s2] // Chop. Commented Sep 14, 2020 at 1:47
• +1. NDSolve and/or ParametricNDSolve are a natural approach to such problems. Commented Sep 14, 2020 at 11:26