# Help me to use DSolve

I am trying to solve the following three coupled ODEs with six boundary conditions:

deq = {pr[x] + a y'[x] + b pr'[x] - c pr''[x] == 0,
e pt[x] + e2 pr'[x] + l pt'[x] + b2 pr[x]/x - dd pt''[x] == 0,
v0 pr'[x] + v0 pr[x]/x - dd y'[x]/x - dd y''[x] == 0, y[L] == y0,
pt[L] == pt0, pr[L] == pr0, y'[L] == v0 pr[L]/dd,
pr'[L] == (pr[L] + v0 y[L])/(2 dd),
pt'[L] == v0 (b2 pr[L] + e pt[L])/(2 dd)} /.
Thread[Rule[{a, b, c, e, e2, l, b2, dd, L, pt0, pr0, v0}, ConstantArray[1, 12]]];

DSolve[deq, {pr, pt, y}, x]


but there is no meaningful output. Could anyone help me to solve the equations? What is wrong in my code?

I don't get solutions to the equations even when I put constants equal to 1:

L = 5;
deq = {pr[x] +  y'[x] + pr'[x] -  pr''[x] == 0,
pt[x] + pr'[x] + l pt'[x] + pr[x]/x - pt''[x] == 0,
pr'[x] + pr[x]/x - y'[x]/x - y''[x] == 0, y[L] == y0,
pt[L] == pt0, pr[L] == pr0, y'[L] == pr[L],
pr'[L] == (pr[L] + y[L]), pt'[L] == pr[L] + pt[L]} /.
Thread[Rule[{a, b, c, e, e2, l, b2, dd, L, pt0, pr0, v0}, ConstantArray[1, 12]]];

DSolve[deq, {pr, pt, y}, x]


I tried the first answer as below:

L = 5;
DI = 1;
f = 1;
g = 1;
ff = 1;
fff = 1;
fo = 1;
gg = 1;
DD = 1;
Y0 = 1;
pt0 = Cos[Pi/6];
pr0 = Sin[Pi/6];
deq = {pr[x] + g y'[x] + gg pr'[x] - DD pr''[x] == 0,
pt[x] + f pr'[x] + ff pt'[x] + fff pr[x]/x - fo pt''[x] == 0,
v0 pr'[x] + v0 pr[x]/x - DI y'[x]/x - DI y''[x] == 0, y[L] == Y0,
pt[L] == pt0, pr[L] == pr0, y'[L] == v0 pr[L]/DI,
pr'[L] == (pr[L] + y[L]), pt'[L] == pr[L] + pt[L]} /.
Thread[Rule[{a, b, c, e, e2, l, b2, dd, L, pt0, pr0, v0, y0},
ConstantArray[1, 13]]];

sol = NDSolve[deq, {pr, pt, y}, {x, 10^-6, 4}];

Plot[Evaluate[{pr[x], pt[x], y[x]} /. sol], {x, 0, 5}, PlotLegends -> Placed[{"pr[x]", "pt[x]", "y[x]"}, {.75, .3}]]


but the output graph shows that pt and pr are not equal to boundary condition at x=l.

• If a symbolic solver returns unevaluated, then Mathematica can't solve the problem. I believe anyone asking for help with such a case should demonstrate that it is reasonable to expect that the problem has a closed form solution. Mathematica is not a magic box that'll spit out a solution to any problem. – Szabolcs Jan 17 '18 at 9:42
• @Szabolcs So true... – Henrik Schumacher Jan 17 '18 at 9:50

Solution by Maple 2017.3 (By magic box)

I chose these to solve:

deq = {pr[x] + a y'[x] + b pr'[x] - c pr''[x] == 0,
e pt[x] + e2 pr'[x] + l pt'[x] + b2 pr[x]/x - dd pt''[x] == 0,
v0 pr'[x] + v0 pr[x]/x - dd y'[x]/x - dd y''[x] == 0, y[L] == y0,
pt[L] == pt0, pr[L] == pr0, y'[L] == v0 pr[L]/dd,
pr'[L] == (pr[L] + v0 y[L])/(2 dd),
pt'[L] == v0 (b2 pr[L] + e pt[L])/(2 dd)} /.
Thread[Rule[{a, b, c, e, e2, l, b2, dd, L, pt0, pr0, v0}, ConstantArray[1,12]]];


Translated from Maple to MMA:

  {pr[x] = (-(1/6))*E^(1 - x)*(-3 + y0) + (1/6)*E^(-2 + 2*x)*(3 + y0),
pt[x] = -((E^(-(5/2) - Sqrt[5]/2 + (1/2)*(1 + Sqrt[5])*x)*(-5*E^2 -
5*Sqrt[5]*E^2))/(10*Sqrt[5])) +
(E^(-(5/2) + Sqrt[5]/2 - (1/2)*(-1 + Sqrt[5])*x)*(-5*E^2 +
5*Sqrt[5]*E^2))/(10*Sqrt[5]) +
(E^(Sqrt[5]*x - (1/2)*(-1 + Sqrt[5])*x)*Integrate[(-(1/6))*E*(-3 + y0)*(-
E^(-((3*z1)/2) - (Sqrt[5]*z1)/2) + E^(-((3*z1)/2) - (Sqrt[5]*z1)/2)/z1) +
((3 + y0)*(2*E^((3*z1)/2 - (Sqrt[5]*z1)/2) + E^((3*z1)/2 -
(Sqrt[5]*z1)/2)/z1))/(6*E^2), {z1, 1, x}])/Sqrt[5] -
Integrate[(-(1/6))*E*(-3 + y0)*(-E^(-((3*z1)/2) + (Sqrt[5]*z1)/2) + E^(-
((3*z1)/2) + (Sqrt[5]*z1)/2)/z1) +
((3 + y0)*(2*E^((3*z1)/2 + (Sqrt[5]*z1)/2) + E^((3*z1)/2 +
(Sqrt[5]*z1)/2)/z1))/(6*E^2), {z1, 1, x}]/(E^((1/2)*(-1 +
Sqrt[5])*x)*Sqrt[5]),
y[x] == y0 + Integrate[-(((1/6)*E^(1 - z1)*(-3 + y0)*z1 - (1/6)*E^(-2 +
2*z1)*(3 + y0)*z1)/z1), {z1, 1, x}]}


And Integrated for x>1:

 solA = {pr[x] = -(1/6) E^(1 - x) (-3 + y0) + 1/6 E^(-2 + 2 x) (3 + y0),
pt[x] = -((
E^(-(5/2) - Sqrt[5]/2 +
1/2 (1 + Sqrt[5]) x) (-5 E^2 - 5 Sqrt[5] E^2))/(10 Sqrt[5])) + (
E^(-(5/2) + Sqrt[5]/2 -
1/2 (-1 + Sqrt[5]) x) (-5 E^2 + 5 Sqrt[5] E^2))/(10 Sqrt[5]) + (
1/(6 Sqrt[5]))
E^(-2 + Sqrt[5] x -
1/2 (-1 + Sqrt[5]) x) ((
2 E^(3 -
1/2 (3 + Sqrt[5]) x) (-1 + E^(
1/2 (3 + Sqrt[5]) (-1 + x))) (-3 + y0))/(3 + Sqrt[5]) - (
4 E^(-(Sqrt[5]/
2)) (-E^(3/2) + E^(1/2 (Sqrt[5] - (-3 + Sqrt[5]) x))) (3 +
y0))/(-3 + Sqrt[
5]) + (3 + y0) (-ExpIntegralEi[1/2 (3 - Sqrt[5])] +
ExpIntegralEi[-(1/2) (-3 + Sqrt[5]) x]) +
E^3 (-3 + y0) (ExpIntegralEi[1/2 (-3 - Sqrt[5])] -
ExpIntegralEi[-(1/2) (3 + Sqrt[5]) x])) -
1/(6 Sqrt[5])
E^(-2 - 1/2 (-1 + Sqrt[5]) x) ((
2 E^(3/2) (-E^((Sqrt[5]/2)) + E^(
1/2 (3 + (-3 + Sqrt[5]) x))) (-3 + y0))/(-3 + Sqrt[5]) + (
4 (-E^(1/2 (3 + Sqrt[5])) + E^(1/2 (3 + Sqrt[5]) x)) (3 + y0))/(
3 + Sqrt[5]) +
E^3 (-3 + y0) (ExpIntegralEi[1/2 (-3 + Sqrt[5])] -
ExpIntegralEi[1/2 (-3 + Sqrt[5]) x]) + (3 +
y0) (-ExpIntegralEi[1/2 (3 + Sqrt[5])] +
ExpIntegralEi[1/2 (3 + Sqrt[5]) x])),
y[x] = y0 +
1/12 (3 + 2 E^(1 - x) (-3 + y0) - 3 y0 + E^(-2 + 2 x) (3 + y0))}


Numerical check:

(* for analitic solution: *)

y0 = 1;
solA /. x -> 4 // N
(* {268.969, 503.393, 135.46} *)

(* for numeric solution:*)

L = 1;
sol = NDSolve[deq, {pr, pt, y}, {x, 1, 4}];
{pr[x], pt[x], y[x]} /. sol /. x -> 4 // N
(*{{268.969, 503.393, 135.46}}*)

• Thank you for your great answer :). Do you think below equation can also be solved by the same method? – Wat Watson Jan 17 '18 at 13:13
• EQ = {y[x] + x y'[x] == v0 pr[x]/Dd + x v0 pr'[x]/Dd, (4 v0 pr0 x + Dd) pr[x] + v0 x^2 y[x]/2 + (2 v0 x^2 pr0 - Dd x) pr'[x] - x^2 Dd pr''[x] == 0, y[L] == 1, pr[L] == pr0, y'[L] == v0 pr[L]/Dd, pr'[L] == v0 (pr0 pr[L] + y[L]/2)/Dd} – Wat Watson Jan 17 '18 at 13:13
• Something seems to be wrong… just try solrule = Unevaluated@(*The last code block*)/. HoldPattern@(a_[x] = b_) :> (a -> Function @@ {x, b}); Simplify[deq /. solrule, x > 1] The 1st equation isn't satisfied. – xzczd Jan 17 '18 at 13:21
• @WatWatson.Maple says an Error:Error, (in dsolve) found differentiated functions with same name but depending on different arguments in the given DE system: {pr(L), pr(x), y(L), y(x)}  – Mariusz Iwaniuk Jan 17 '18 at 13:24
• @xzczd.I edited my answer. Maple's simplify command something does not work properly or converting from MathML is weird. – Mariusz Iwaniuk Jan 17 '18 at 13:53
L = 5;
deq = {pr[x] + y'[x] + pr'[x] - pr''[x] == 0,
pt[x] + pr'[x] + l pt'[x] + pr[x]/x - pt''[x] == 0,
pr'[x] + pr[x]/x - y'[x]/x - y''[x] == 0, y[L] == y0, pt[L] == pt0,
pr[L] == pr0 , y'[L] == pr[L], pr'[L] == (pr[L] + y[L]),
pt'[L] == pr[L] + pt[L] } /.
Thread[Rule[{a, b, c, e, e2, l, b2, dd, L, pt0, pr0, v0, y0},
ConstantArray[1, 13]]];

sol = NDSolve[deq, {pr, pt, y}, {x, 10^-6, 4}];

LogPlot[Evaluate[{pr[x], pt[x], y[x]} /. sol], {x, 10^-6, 4},
PlotLegends -> Placed[{"pr[x]", "pt[x]", "y[x]"}, {.75, .3}]]


• Thank you for your answer. But I prefer to solve them without ndsolve to understand the effect of each parameter. Do you think it could be possible to do that or not? – Wat Watson Jan 17 '18 at 6:11
• Actually, I need to find a functionality for the solutions. – Wat Watson Jan 17 '18 at 6:13
• Since DSolve cannot solve the system with numeric values for all parameters, it is extremely unlikely that a solution can be found for symbolic parameters. – Bob Hanlon Jan 17 '18 at 6:20
• What do you mean by "numeric values for all parameters"? Can DSolve find the functionality if I determine some parameters and not all of them? – Wat Watson Jan 17 '18 at 6:25
• The Thread[Rule[{a, b, c, e, e2, l, b2, dd, L, pt0, pr0, v0, y0}, ConstantArray[1, 13]] sets all of the parameters to 1, i.e., assigns "numeric values for all parameters". Even with these assignments, DSolve cannot solve the simplified system. Complicating the problem by leaving any of the parameters as symbolic is not going to help DSolve. – Bob Hanlon Jan 17 '18 at 6:33

If you try to find the general solution of your equations without any boundary conditions

DSolve[{pr[x] + Derivative[1][pr][x] +Derivative[1][y][x] - (pr^\[Prime]\[Prime])[x] == 0,
pr[x]/x + pt[x] + Derivative[1][pr][x] +Derivative[1][pt][x] - (pt^\[Prime]\[Prime])[x] == 0,
pr[x]/x + Derivative[1][pr][x] - Derivative[1][y][x]/x - (y^\[Prime]\[Prime])[x] == 0 }
, {pr, pt, y}, x]

(*DSolve[{pr[x] + Derivative[1][pr][x] +Derivative[1][y][x] - (pr^\[Prime]\[Prime])[x] == 0,
pr[x]/x + pt[x] + Derivative[1][pr][x] +Derivative[1][pt][x] - (pt^\[Prime]\[Prime])[x] == 0,
pr[x]/x + Derivative[1][pr][x] - Derivative[1][y][x]/x - (y^\[Prime]\[Prime])[x] == 0}, {pr, pt, y}, x]*)


MMA can't solve it! So any further activity using DSolve for the complete system seems to be useless!!!

Let's try some simplication: The first and the last ode only depends on y[x] , pr[x] and can be separated:

Derivative[1][y][x] = -pr[x] -Derivative[1][pr][x] + (pr^\[Prime]\[Prime])[x] (* first ode *)


substituting in the third ode

oder=2 pr[x] + (1 + 2 x) Derivative[1][pr][x] - (pr^\[Prime]\[Prime])[x] +x(pr^\Prime]\[Prime])[x] - x\!$$\*SuperscriptBox[\(pr$$,TagBox[RowBox[{"(","3", ")"}],Derivative],MultilineFunction->None]\)[x]


gives ode in pr which can be solved

ergpr = DSolve[oder == 0, pr, x][[1]]
(* {pr -> Function[{x},E^-x C[1] + E^(2 x) C[2] +3 E^-x C[3] (1/3 E^(3 x) ExpIntegralEi[-2 x] - ExpIntegralEi[x]/3)]} *)


With this result y'[x]is known and y[x] can be evaluated

ergy= DSolve[y'[x] == -pr[x] -Derivative[1][pr][x] + (pr^\[Prime]\[Prime])[x] /.ergpr, y,x][[1]]
(*{y -> Function[{x}, -E^-x C[1] + 1/2 E^(2 x) C[2] + C[4] +
1/2 E^(2 x) C[3] ExpIntegralEi[-2 x] +
E^-x C[3] ExpIntegralEi[x] - 1/2 C[3] Log[-2 x] + 2 C[3] Log[x]]}*)


Unfortunately the remaining second ode

odept = -E^-x C[1] + 2 E^(2 x) C[2] +3 E^(2 x) C[3] ExpIntegralEi[-2 x] +(E^-x C[1] + E^(2 x) C[2] +3 E^-x C[3] (1/3 E^(3 x) ExpIntegralEi[-2 x] -ExpIntegralEi[x]/3))/x -3 E^-x C[3] (1/3 E^(3 x) ExpIntegralEi[-2 x] -ExpIntegralEi[x]/3) +pt[x] + Derivative[1][pt][x] - (pt^\[Prime]\[Prime])[x] == 0


can not be solved analytically by MMA!

Adaption of the boundary condition y, pr perhaps simplifies the last equation...