The DE has problems at y[x] == 0
as noted, but also at x == 1
, where there is a pole of order 2. This suggests that numerical issues near the singularities could cause substantial error to accumulate. One can use Piecewise
to substitute the limiting value at a discontinuity, but this won't cure any numerical instability in the neighborhood of the discontinuity.
One can remove the 1/0
errors with "EquationSimplification" -> "Residual"
, but the numerical trouble persists. The code below results in a plot similar to, if not worse than, the OP's workaround.
sol = NDSolve[{DE == 0, y[0] == 0, y'[0] == -4, y''[0] == 24},
y, {x, 0, 10}, Method -> {"EquationSimplification" -> "Residual"}]
Plot[{(-4 x)/(1 + 3 x), y[x] /. First[sol]}, {x, 0, 10},
PlotStyle -> {AbsoluteThickness[5], Automatic}]

One can see that the second derivative is not calculated very well near x == 0
by this method:
foo = D[(-4 x)/(1 + 3 x), x, x];
Plot[{y''[x] /. First[sol], foo}, {x, 0, 0.001}, PlotRange -> All]

Using Piecewise
helps near x == 0
, but it fails at x == 1
:
sys = {y''[x] ==
Piecewise[{{y''[x] /. First@Solve[DE == 0, y''[x]] // Expand //
Simplify, y[x] != 0}}, 24], y[0] == 0, y'[0] == -4};
sol = NDSolve[sys, y, {x, 0, 10}];
NDSolve::ndsz: At x == 1.00008646582443`, step size is effectively zero; singularity or stiff system suspected. >>
Plot[{(-4 x)/(1 + 3 x), y[x] /. First[sol]}, {x, 0, 1.},
PlotStyle -> {AbsoluteThickness[5], Automatic}]

To get a solution that crosses the discontinuity is difficult. The numerical evidence* suggests the discontinuity is unstable and virtually impossible without some preconditioning. The best thing I have to offer is to start the integration at x == 1
. While there are some numerical issues near x == 0
, integrating toward it seems stable.
(*I tried WhenEvent[x == 1, "CrossDiscontinuity"]
, increasing WorkingPrecision
, and such stuff. Whenever it integrated past x == 1
, there was a large error.)
To start the integration at x == 1
, I tried to analyze the DE and come up with a solution (pretending not to know the OP's desired outcome). It turns out there are two choices for y[1]
, so that didn't work. And for the desired y[1] == -1
, the choice for y'[1]
appears to be free, so again I had to use the OP's formula to determine the initial value. I determined these constraints by examining the series expansion of the DE. The default value for y''[1]
is also determined this way.
SeriesCoefficient[DE, {x, 1, -2}] == 0 // Solve
SeriesCoefficient[DE, {x, 1, -1}] == 0 /. y[1] -> -1 // Solve
SeriesCoefficient[DE, {x, 1, 0}] == 0 /. y[1] -> -1 // Solve
(*
{{y[1] -> -2}, {y[1] -> -2}, {y[1] -> -1}}
{{}}
{{y''[1] -> y'[1] (-1 + 2 y'[1])}}
*)
Note that "StiffnessSwitching"
is needed to step through the muck around x == 1
:
ypp[x0_?NumericQ, y0_?NumericQ, yp0_?NumericQ] = Piecewise[
{{y''[x] /. First@Solve[DE == 0, y''[x]] /. {x -> x0, y[x] -> y0,
y'[x] -> yp0}, x0 != 1}},
yp0 (-1 + 2 yp0)];
sol = NDSolve[{y''[x] == ypp[x, y[x], y'[x]], y[1] == -1,
y'[1] == -1/4}, y, {x, 0, 10}, Method -> "StiffnessSwitching"]
Plot[{(-4 x)/(1 + 3 x), y[x] /. First[sol]}, {x, 0, 10},
PlotStyle -> {AbsoluteThickness[5], Automatic}]

y=-2
. Another isy=(-4*((x^4-2*x^2+1)/(x^4+8*x^3+18*x^2+8*x+1))^(1/2)*x^4-20*((x^4-2*x^2+1)/(x^4+8*x^3+18*x^2+8*x+1))^(1/2)*x^3-4*x^4-24*((x^4-2*x^2+1)/(x^4+8*x^3+18*x^2+8*x+1))^(1/2)*x^2-16*x^3-20*((x^4-2*x^2+1)/(x^4+8*x^3+18*x^2+8*x+1))^(1/2)*x-32*x^2-4*((x^4-2*x^2+1)/(x^4+8*x^3+18*x^2+8*x+1))^(1/2)-16*x-4)/(3*x^4+8*x^3+14*x^2+8*x+3)
as well as the.. $\endgroup$ – Nasser Sep 23 '14 at 8:54y=-4*x/(1+3*x)
. The point is, these are singular solutions. There is a singularity in your ode, which is why NDSolve gives a 1/0 error. So, there is nothing wrong withNDSolve
. Notice thatDSolve
did not solve this. $\endgroup$ – Nasser Sep 23 '14 at 8:55