A single call to DSolve
can, in fact, return the solution to this pair of ODEs, but it is not pretty. For convenience, define
eq = {y0[x] + y0''[x] == 0, y1[x] + y1''[x] == x y0[x]};
Then, as noted in the question, two calls to DSolve
yield
ry0 = DSolve[First@eq, {y0[x]}, x] // Flatten
ry1 = DSolve[Last@eq /. %, y1[x], x] // Flatten // Simplify
(* {y0[x] -> C[1] Cos[x] + C[2] Sin[x]} *)
(* {y1[x] -> 1/8 ((2 x C[1] + C[2] - 2 x^2 C[2] + 8 C[3]) Cos[x] +
((-1 + 2 x^2) C[1] + 2 x C[2] + 8 C[4]) Sin[x])} *)
To obtain the equivalent result with one call to DSolve
, combine the two equations into one
syt = Last@SolveValues[Last@eq, y0[x]]
First@eq /. y0 -> Function[{x}, Evaluate@%];
eq4 = Simplify[x^3 #] & /@ %
(* (y1[x] + y1''[x]/x *)
(* (2 + x^2)*y1[x] - 2*x*y1'[x] + 2*y1''[x] + 2*x^2*y1''[x] - 2*x*y1'''[x] + x^2*y1''''[x] == 0 *)
and solve eq4
for y1
.
DSolveValue[eq4, y1[x], x, GeneratedParameters -> CC];
cf[e_] := LeafCount[e] + 100 Count[e, _ArcTan, {0, Infinity}] +
100 Count[e, _Sqrt, {0, Infinity}] + 100 Count[e, _Rational, {0, Infinity}];
FullSimplify[%%, x > 1, ComplexityFunction -> cf];
sy1 = Simplify[% // ExpToTrig // TrigExpand]
(* (CC[1] + CC[2] + (-1 + x) (-I ((1 + I) + x) CC[3] + ((1 + I) + I x) CC[4])) Cos[x] -
I (CC[1] - CC[2] + I (-1 + x) (((1 + I) + x) CC[3] + ((1 - I) + x) CC[4])) Sin[x] *)
To obtain the solution for y0
, back-substitute.
sy0 = Simplify[syt /. y1 -> Function[{x}, Evaluate@sy1]]
(* 4 (CC[3] + CC[4]) Cos[x] + 4 I (CC[3] - CC[4]) Sin[x] *)
Finally, the equivalence of this solution to that obtained with two calls to DSolve
can be seen from
Simplify[{sy0 - (y0[x] /. ry0), sy1 - (y1[x] /. ry1)}];
Last@CoefficientArrays[%, {Cos[x], Sin[x]}];
Solve[Thread[Flatten@% == 0], Array[CC, 4]] // Flatten
(* {CC[1] -> 1/16 ((2 + I) C[1] - (1 - 2 I) C[2] + 8 C[3] + 8 I C[4]),
CC[2] -> -(1/16) I ((1 + 2 I) C[1] + (2 - I) C[2] + 8 I C[3] + 8 C[4]),
CC[3] -> 1/8 (C[1] - I C[2]),
CC[4] -> 1/8 (C[1] + I C[2])} *)
MatrixExp[]
$\endgroup$