# DSOLVE for tenth order boundary value problem

I am currently using Mathematica 9 to find solutions for higher order ordinary differential equations. The following was my input for a tenth order boundary value problem (BVP):

DSolve[{y''''''''''[x] == (E^-x)*(y[x])^2, y == 1, y'' == 1,
y'''' == 1, y'''''' == 1, y'''''''' == 1, y == E,
y'' == E, y'''' == E, y'''''' == E, y'''''''' == E},
y[x], x]


DSolve returned unevaluated. So, I tried for a numerical solution with slight modification by giving the range for x:

NDSolve[{y''''''''''[x] == (E^-x)*(y[x])^2, y == 1, y'' == 1,
y'''' == 1, y'''''' == 1, y'''''''' == 1, y == E,
y'' == E, y'''' == E, y'''''' == E, y'''''''' == E},
y[x], {x, 0, 1}]


But, it gave me the result in terms of an interpolating function, so fine until here. The third time I used the last input but this time with DSolve that is :

DSolve[{y''''''''''[x] == (E^-x)*(y[x])^2, y == 1, y'' == 1,
y'''' == 1, y'''''' == 1, y'''''''' == 1, y == E,
y'' == E, y'''' == E, y'''''' == E, y'''''''' == E},
y[x], {x,0,1}]


However, I get the error

DSolve::dsvar: 0 cannot be used as a variable. and the output was same as input.


I know the exact solution for the BVP which is $y(x)=e^x$, But, how do I obtain it using Mathematica? What am I doing wrong in the first and the third case when I am using DSolve for getting exact solution instead of a numeric one. Is there some error or there is some limitation in the software?

Again let me clarify, i.e. it is possible to find a numerical solution using NDSolve[] for this example. So it's true that equation has an approximate solution which can be given in form of some interpolating function. You can check this by simply putting in my NDSolve input in Mathematica. The point is, this equation has a solution which is $y(x)=e^x$. Here is a link to the research paper where this is mentioned (please check example 3.1. Note that in the paper, there is a small mistake: boundary condition should be $y(8)(1)=e$ instead of $y(8)(0)=e$ since it has already been mentioned earlier that $y(8)(0)=1$. So there cannot be two different values at one point. Rest of it is fine and one can see that it is clearly written that the exact solution is $y(x)=e^x$.).

My question is simple

1. If there is an exact solution to the problem above (which I've already told that there is an exact analytic solution) then how to get it through Mathematica.
2. Or there is some limitation of the software due to which we can have numerical solutions to say, problems of tenth or higher order BVPs and no analytic solutions.

Please, be specific and, yes, I have a fair understanding that if an equation has a general solution that same one can be evaluated given boundary conditions.

• Try NDSOlve it works: !Mathematica graphics (this is non-linear ODE. Very hard to find analytical solution for) and try to use the D[y[x],{x,n}] form when you have such large ODE. Easier to read than y'''''''''''''[x] and you can't use {x,0,1} with DSolve, only with NDSolve Aug 28 '13 at 14:29
• @Nasser he tried NDSolve, and it worked. His question is why doesn't DSolve work. On your other points, I agree, Mathematica is not likely to be able to do this automatically. Aug 28 '13 at 14:34
• you will notice that in case of NDSolve default range is {x,0,1} if you need more only than u shall specify. Aug 29 '13 at 18:57

I do not think the problem is with the derivative order here. Mathematica's DSolve first solves general problem and then looks for the constants. So there is no harm here to check what kind of complexity we dealing here with at very first steps:

1st order - solvable

DSolve[D[y[x], {x, 1}] == (E^-x)*(y[x])^2, y[x], x]


{{y[x] -> -(E^x/(-1 + E^x C1))}}

2nd order - unsolvable

DSolve[D[y[x], {x, 2}] == (E^-x)*(y[x])^2, y[x], x]


Further derivative order increase will just increase complexity of this solution. And just as proof of concept - here is 10th order general case of solvable thing:

DSolve[D[y[x], {x, 10}] == x + y[x], y[x], x] // TraditionalForm And here is how you do just 5 initial conditions substitution, no need to type '''''... etc.

EQ = {D[y[x], {x, 10}] == x + y[x],
Table[D[y[x], {x, n}] == 0, {n, 0, 4, 1}] /. x -> 0}


Just with 5 you get a million page solution and if you would do all 9 then - brace your self. I am trying to say the problem here is not in the order of the derivative, but in underlying complexity of your equation already at order 2. Not all deferential equations are intagrable by standard methods or at all. Do not blame the order. While equation admits simple particular solution, the ODE itself is non-liner. The non-linear ODE do not have universal solving techniques, this maybe the reason it is hard to find that Exp solution.

DSolve[EQ, y[x], x] // TraditionalForm I tried to investigate a bit and I am getting a little contradictory result that I will explain as under. You have specified that the result you obtain is $E^x$ but on following the below steps, I am seeing some variation.

ifun = First[
y /. NDSolve[{y''''''''''[x] == (E^-x)*(y[x])^2, y == 1,
y'' == 1, y'''' == 1, y'''''' == 1, y'''''''' == 1,
y == E, y'' == E, y'''' == E, y'''''' == E,
y'''''''' == E}, y, x]]


InterpolatingFunction[{{0.,1.}},<>]

GraphicsRow[{Plot[ifun[x], {x, -1., 1.}], Plot[ifun[x], {x, 0., 5.}]}]


I tried if I can find some functionality that one can determine equation for some data and I found this Lagrange_polynomial and if you are really interested in finding out, you can program it. In fact it must be present n M already but I am not aware of it. Now I tried to find some interpolating polynomial just to try if I will get some equation like this.

ip = InterpolatingPolynomial[Table[ifun[x], {x, 10}], x]
sip = Simplify[ip]


-3.93862*10^37 + 1.11928*10^38 x - 1.28662*10^38 x^2 + 8.00427*10^37 x^3 - 3.00906*10^37 x^4 + 7.15674*10^36 x^5 - 1.08482*10^36 x^6 + 1.01649*10^35 x^7 - 5.36862*10^33 x^8 + 1.22275*10^32 x^9 On plotting this,

Plot[sip, {x, -263.437, 263.437}]


we get the same graph, which shows that it was a valid step. Now why I am saying this is because graph of $y=E^x$ never goes in -y. As Limit[E^x, x -> -Infinity] return 0. The positive rise of your equation rises like $E^x$ but not so in case of -y. I believe I might be making mistake(as you have published result), but I am not seeing where. Why results are not matching your conclusion.?

• ifun is only defined from 0 to 1 by NDSolve. Sep 24 '16 at 14:47