# DSolve inconsistently able to evaluate expression

I am trying to use DSolve to evaluate a coupled system of ODEs as follows

DSolve[{Y0[X] + Y0''[X] == 0, Y1[X] + Y1''[X] == X Y0[X]}, {Y0[X], Y1[X]}, X]


However, when I do so DSolve refuses to evaluate and returns nothing. When I instead evaluate the equations one at a time, as

sol = DSolve[{Y0[X] + Y0''[X] == 0}, {Y0[X]}, X]
DSolve[{Y1[X] + Y1''[X] == X Y0[X]} /. sol, Y1[X], X]


I am given the correct analytic forms of the solution. This behavior holds for evaluating with or without boundary conditions. A set of boundary conditions that yields analytical results would be $$y_0(0)=1$$ and $$y_0'(0) = y_1(0) = y_1'(0) = 0$$, which evaluate following the same behavior as above. Does anyone know why this behavior exists? In this example it is fairly straight forward to do each step sequentially but as I plan to do something similar with much more complicated expressions I'd like to work the edge cases out ahead of time.

• There was a question like this in the past few years. (I remember answering or investigating one, at least.). When I get a chance I’ll look for it, or maybe someone else will find it first. I don’t think the other one was linear, though. Sep 14 at 20:34
• As a workaround transform your ode to a first order system and solve it using MatrixExp[] Sep 15 at 8:32
• Unfortunately while in this case such a transformation is possible, in the general case that I am hoping to solve next it would not be. Therefore, I was hoping to try and solve this behavior on this simpler problem before proceeding to work on my next problem. Sep 15 at 21:08
• What is your more general system of ODEs? A reader may be able to solve it with a bit of work. Sep 18 at 12:44
• I don't actually have it myself - was trying to provide some assistance to a friend but I was in general looking for a more systematic way to handle these kind of problems. Specifically, handing solving asympototic expansions of differential equations, which is what these two equations represent the first two terms of. Sep 19 at 20:11

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])} *)

• Why is it that you need to go through this process instead of DSolve working as posted in the first part? Sep 19 at 13:05
• @JamesVR It seemed to me that you were asking how to obtain the solution without calling DSolve twice. I provided this answer in the hope that this might help with the more complex case that you have not shared. There are many "simple" systems of differential equations that DSolve cannot treat. You might wish to inform Wolfram, Inc of this case. By the way, I am told that Maple is more powerful in solving differential equations symbolically, but I have no personal experience with it. Sep 19 at 14:38