# Is it possible to fit such a large range plot in mathematica?

I am trying to solve the coupled differential equation numerically with Mathematica. But the range of values are large so mathematica cannot give correct plot. Here is the code:

a = 4.75388*10^26;
b = 5.424*10^-3;
d = 4.75388*10^20;
{X, Y} = {x, y} /.
NDSolve[{x'[
z] == -((a/z) (x[z] - b*z^(3/2) E^(-z)) (BesselK[1, z]/BesselK[2, z])),
y'[z] == ((d/z) (x[z] - b *z^(3/2) E^(-z)) (BesselK[1, z]/
BesselK[2, z]) - (a *z/4) (BesselK[1, z]) y[z]),
x[0.1] == 1.552*10^-4, y[0.1] == 10^(-9)}, {x, y}, {z,0.1,100}] //
FullSimplify // First
LogLogPlot[{X[z], Y[z]}, {z, 0.1, 100}, PlotRange -> All]


I am getting an Interpolating function whose domain is restricted to only {0.1,0.1}. I think this is the main reason of all the errors

• This may be an obvious question: What about normalization of solution (to between 0 and 1) and then plotting it? Does that help? – dearN Sep 17 '19 at 11:48
• When I run your code the first error I get says NDSolve::ndlim: Range specification z is not of the form {x, xend} or {x, xmin, xmax}. – N.J.Evans Sep 17 '19 at 12:23
• Are we sure the issue is the size of the numbers and not the 5 errors that are given on evaluation? – ktm Sep 17 '19 at 12:57
• When I run the code, NDSolve fails to take the first step (integration step); error NDSolve::ndsz. If it did succeed, the plot would fail because X is not a function, but a replacement Rule; hence X[z], while syntactically allowable, is semantically meaningless. If NDSolve can be fixed, try x[z] /. X as shown in the documentation. – Michael E2 Sep 18 '19 at 11:28

Geneal approach to problems like this is to attempt to regularize your DE. Naive approach in this case would be to replace $$x$$ with $$x1=x/a$$ and $$y$$ with $$y1=y/d$$:

a = 4.75388*10^26;
b = 5.424*10^-3;
d = 4.75388*10^20;

{X, Y} = {x1[#]*a &, y1[#]*d &} /. NDSolve[
{
x1'[z] == -((BesselK[1, z] (-b E^-z z^(3/2) + a x1[z]))/(
z BesselK[2, z]))
, y1'[z] == (BesselK[1, z] (-b E^-z z^(3/2) + a x1[z]))/(
z BesselK[2, z]) - (a z BesselK[1, z] y1[z])/4
, x1[0.1] == 1.552*10^-4/a
, y1[0.1] == 10^(-9)/d
}
, {x1, y1}
, {z, 0.1, 100}
] // First
LogLogPlot[{X[z], Y[z]}, {z, 0.1, 100}, PlotRange -> All]


Edit:

Resulting $$y(z)$$ does satisfy boundary conditions:

In[]:= Y[0.1]

Out[]= 1.*10^-9


But it might be not so obvious from picture above.

Plot[Y[z], {z, 0.1, 0.1 + 10^-10}, PlotRange -> All]


There is also an interval in which $$x(z)$$ becomes negative that is displayed somewhat misleading on the LogLogPlot

• it is a very approach but why not the y[z] initially starting from near the 10^-9 point since it is one of the given boundary conditions to solve this differential equation – user105697 Sep 18 '19 at 10:04
• @user105697 edited reply to your concern into my answer. – Markhaim Sep 18 '19 at 10:36
• Thank you for your effort. I think i have to revisit my input equations carefully. – user105697 Sep 18 '19 at 11:29
• BTW, a more straightforward way to get the functions X, Y is to replace {x1[#]*a &, y1[#]*d &} /. NDSolve[args...] // First with NDSolveValue[args...]. – Michael E2 Sep 18 '19 at 11:48