# How do I find out what a Mathematica function is doing (which numerical method it uses and how)?

Long story short, I am working on a project which entails fixing a bug I have found in Mathematica (8 and 10). In general, I would like to know what a specific Mathematica function is actually doing when it is called. That is, what is the numerical method implemented by a Mathematica function.

For my specific case, I am trying to figure out how Plot() takes a InterpolatingFunction and fits a curve to it. I have found cases where the fit is dependent on the interval which is plotted. If I expand or contract the plotting domain I get erroneous curve fitting in a localized region. It seems that Plot() selects a subset of data points from the InterpolatingFunction and fits the data. This causes error (though rare and easily avoided) for a specific curve I am plotting.

I am not asking for a solution to this problem, I just want to know what the Plot() function does when I give it an InterpolatingFunction. In general, I would like to obtain this information for any function. It seems like Wolfram would document this information somewhere, but nothing this specific can be found in the documentation center.

Here is what I am investigating:

para = {Γ -> (0.029), g0 -> (9.6*10^-5),
ntr -> (4*10^17), β -> (1.8*10^-4), r -> (3.1*10^-9),
B -> 1/(nt r), nt -> 1/(Γ ph g0) + ntr,
ph -> (1.25*10^-12), Jt -> nt/r, st -> Γ ph nt/r,
pj -> 1, ppj -> 1};

JC[t_] := 0.9 Jt + 0.3 Jt UnitStep[t - 20*10^-9] UnitStep[Cos[t 7*10^8]] //.     para;

eqs = {n'[t] == JC[t] - g0 (n[t] - ntr) S[t] - B (n[t])^2,
S'[t] == Γ g0 (n[t] - ntr) S[t] - S[t]/ ph + Γ β B (n[t])^2, n == 0,
S == 0} //. para;

sol = NDSolve[eqs, {n, S}, {t, 0, 60*10^-9}, MaxSteps -> 10^6]

p1 = Plot[Evaluate[S[t]*10^-12 /. sol], {t, 20 10^-9, 50 10^-9}]


You will see three pulses of decaying oscillations. All should be approximately identical, but the third have a linear break in oscillations. You can prevent this by changing the right endpoint of plotting domain: 50 10^-9 -> 49 10^-9 or 51 10^-9.

I have extracted the data points and found them to provide adequate information to fit an accurate curve. Conceptually, I know why there is an issue, but I need to know precisely what Mathematica has decided to do when fitting the curve for a given interval.

Again, I am not looking for someone to do my project for me, just need to know where I can find the information.

Note: setting PerformanceGoal -> "Quality" does nothing.

• I'm sorry, but if it's not documented, it's not documented; you can either try asking WRI (and/or the user community), or reverse-engineering the behavior yourself. But there is no general answer to this, which is a frustration for many of us, as the documentation is often less than we would hope for. For your specific problem there is a known undocumented option that will be highly relevant: setting the refinement control value to 1 degree solves the problem you observe. – Oleksandr R. Nov 30 '14 at 1:21
• @OleksandrR.Thank your for this, it is very relevant. – Sean Bearden Nov 30 '14 at 1:25

With the functions given in the question above:

Pressing the plus sign in the interpolation function object will expand the description to include the method.

sol = NDSolve[eqs, {n, S}, {t, 0, 60*10^-9}, MaxSteps -> 10^6] The problem is not the iterator range but rather the initial number of PlotPoints. PlotPoints -> n specifies the total number of initial sample points to use. If PlotPoints is too small or just starts with a poorly performing set of initial points then some features can be missed. Adaptive procedures controlled by MaxRecursion are used to choose more sample points.

Plot[Evaluate[S[t]*10^-12 /. sol],
{t, 20 10^-9, 50 10^-9},
ImageSize -> 300,
PlotLabel ->
StringForm["PlotPoints -> ", #],
PlotPoints -> #] & /@
{Automatic, 25, 50, 51} // Column • Thank you very much! This gives me a lot of information with which to work. – Sean Bearden Nov 29 '14 at 23:48
• I have noticed that I cannot reproduce the plot by lowering PlotPoints. I assume this is because of MaxRecursion or some other option. How do I find the value to which an option has automatically been set? – Sean Bearden Nov 29 '14 at 23:49
• @Sean I think you should ask separate question regarding reproducing the plot. One question should be generally related only to one problem. Otherwise one could ask additional questions and subquestions forever, and that's not fair to answerers. :) – Adrian Nov 30 '14 at 0:08
• @Adrian Sorry about that! Thanks for the advice. – Sean Bearden Nov 30 '14 at 0:10
• The page on Options for Graphics reference.wolfram.com/language/tutorial/Options.html states that the default for PlotPoints is 50. However, since 25 worked better than Automatic for your plot, this seems unusual. In any event, if your plot is exhibiting unexpected behavior generally you would want to increase rather than decrease PlotPoints. When in doubt, increase the number of PlotPoints to see if there is a change. The cost is that increased PlotPoints slows down the plotting, particularly with computationally intensive functions. – Bob Hanlon Nov 30 '14 at 0:14