# NDSolve for coupled differential equations

I am trying to solve two first order coupled differential equations for the variables yD and yN given as,

sol1 =
NDSolve[
{{yD'[x] == -((1.81*10^11 (yD[x]^2 - yN[x]^2))/x^2),
yN'[x] == -2.35 *10^9 x (yN[x] - 0.0094 E^(-3 x/4) x^(3/2)) +
(1.81*10^11 (yD[x]^2 - yN[x]^2))/x^2 -
(1.81*10^7 (yN[x]^2 - (0.0094 E^(-3 x/4) x^(3/2))^2))/x^2},
{yD[1] == 0.00107, yN[1] == 0.00445}},
{yD, yN}, {x, 1, 30},
Method -> {StiffnessSwitching,Method -> {ExplicitRungeKutta, Automatic}},
WorkingPrecision -> 50,
MaxSteps -> 50000,
AccuracyGoal -> 10,
PrecisionGoal -> 10]


I get an error message

The precision of the differential equation is less than WorkingPrecision (50.).

Increasing the WorkingPrecision further or using InterpolationOrder -> All does not help.

The plot for the solutions

LogLogPlot[
{Evaluate[yD[x] /. sol1], Evaluate[yN[x] /. sol1],}, {x, 1,30},
PlotRange -> {{1, 1000}, {10^-14, 1}}, PlotStyle -> {Gray, Orange}]


shows very wiggling behavior for the yN[x] solution, however the yD[x] solution is very smooth. Can anyone help me in getting rid of the error message and finding a stable, smooth solution for yN[x] also.

• Your use of constants like 1.81*10^11 and 0.0094 run counter against WorkingPrecision -> 50; if you must, run Rationalize[] or SetPrecision[] over these inexact numbers. Commented Jan 31, 2017 at 8:32
• Thank you for your suggestion. Rationalize[ ] worked fine and the error message is gone now. But I still have the same wiggling behavior for yN[x] solution. Is there a problem in numerical interpolation of the solution. Commented Jan 31, 2017 at 14:36

You seemed to have a few issues in your code.

1. Your numbers were expressed in machine precision, which prevented the use of higher precision even when WorkingPrecision was set. This was pointed out in comments as well.
2. Furthermore, your Method directive should have had names of methods as strings, rather than as symbols.

Once those are fixed, and the WorkingPrecision increased appropriately, you can get a smooth plot for both functions. In fact, you get pretty much the same plot:

sol1 = NDSolve[{
yD'[x] == -((181*10^9 (yD[x]^2 - yN[x]^2))/x^2),
yN'[x] == -235*^7 x (yN[x] - 94*^-4 E^(-3 x/4) x^(3/2)) + (181*^9 (yD[x]^2 - yN[x]^2))/
x^2 - (181*^5 (yN[x]^2 - (94*^-4 E^(-3 x/4) x^(3/2))^2))/x^2,
yD[1] == 107*^-5, yN[1] == 445*^-5},
{yD, yN}, {x, 1, 30},
Method -> {"StiffnessSwitching", Method -> {"ExplicitRungeKutta", Automatic}},
WorkingPrecision -> 35
];

LogLogPlot[
Evaluate[{yD[x], yN[x]} /. sol1],
{x, 1, 30},
PlotStyle -> {Gray, Orange}, PlotRange -> All
]