Your title indicates that you want to compute the 2nd derivative of f at 0. Yet your ODEs have no f. Did you mean that you want the 2nd derivative of U at 0?

You must not have transcribed your code correctly to your Question. You say that sub_var is the list that you present. Yet you have the line

for sub_var in subs do

which indicates that subs is the list and sub_var iterates over that list's entries.

@rallezet Yes, that would be fairly easy. If you tell me the coordinates of the center of the ball as a function of t, I'll do it.

@tomleslie Lower residuals are not necessarily better in the real-number setting. Consider the function

f:= x-> ((x+2)^3*(x-3)^2 + eps)*exp(-2*x)

where, as usual, eps is a small positive number. To be concrete, let eps = 10^(1-Digits). It's obvious that f has a unique real root, and that that root is slightly less than -2. Indeed, fsolve returns -2. as the root. And it's obvious that 0 < f(3) < f(-2), so 3 has a lower residual.

@Vaishnavi You need to execute the entire worksheet. Simply executing the fsolve command is not enough. One way to execute the entire worksheet is to press !!! on the toolbar.

@Adam Ledger I kinda agree with you about that. However, long-standing tradition of at least 60 years, in all computer languages that I'm aware of, is that for hexadecimal (and only for hexadecimal) the letters A-F are used for digits over 9. It's important for the human reading of computer memory dumps that the more-compact one-character-per-digit form be used. Hexadecimal is the only base that facilitates that dump reading. That being said, it's essentially only an accident that the BBP algorithm works in base 16. The meta-algorithm that derives the summation formulas might return a formula in any other base. I'm sure that they'd publish a base-10 formula if they could find one.

@vv It is not worth it to increase the precision by increasing the value of Digits. Hardware-float computation is in my experience 20 - 30 times faster than software float (at the same precision). So, if the hardware floats give you only 1 correct digit, that's still better than using software floats to get, say, 9 correct digits. To get the next digit, it's better to increase d by 1 and do the hardware-float computation again. And, we should be able to get more than 1 correct digit with the hardware floats.

A hardware float has a 53-bit mantissa. One bit is the sign bit. That leaves 52-bits of precision, which is 13 hexadecimal digits. So that is the theoretical limit of accuracy of this method. So I don't know why there was talk of d+20 a few messages back.

@Adam Ledger The conversion of the final computed value into hexadecimal characters is clunky, and I hope to change it. The command convert(x, binary) takes a numeric value x---either an integer or a float---and converts it into a base-10 number that only uses the digits 0 and 1, thereby giving a representation of x in base 2. The x is always a float between 0 and 1 in this program. The "%38.36f" converts that float to a string of characters with 36 characters (all 0s and 1s in this case) after the decimal point and 1 character (necessarily 0 in this case) before the decimal point. Those characters plus the decimal point itself add up to 38 characters.

@sand15 I see no practical difference between your formula and the formula that the OP was trying to transcribe, which is

add(f(randvals[i]), i = 1 .. numsamples)*(rhs(rhs(r))-lhs(rhs(r)))/numsamples.

What practical difference do you see between this and your

add(F~(S) * (b-a)/N ) ??

The key thing here being that Tom capitalized Sum.

Would you please post a worksheet showing both your code and the error message? I believe that you may have incorrectly transcribed the code. Or maybe it is an issue of 2D Input being misinterpreted, which happens all too often.

@acer First off, I think that your idea is superb; so I hope that this doesn't sound too critical. I just think that your example is filled with "magic numbers" that make putting it into practice seem easier than it actually is. Indeed, your -2, 2, 0.5, 0.01, 201, and 9 are all magic numbers. (The 9 is a hidden magic number, being the number of visible gridlines in each dimension.) They're magic because irem(201 - 1, 9 - 1) = 0 and (2 - (-2))/(9 - 1) = 0.5 and frem(-2, 0.5) = 0 and (2 - (-2))/(201 - 1) > 0.01.

Let's suppose that the original plot is

plot3d(F(x,y), x= a..b, y= c..d, grid= [N,M], style= surface, ...),

and that we want n by m of those gridlines to be visible. I impose the further reasonable conditions (as you did implicitly) that we want the gridlines evenly spaced in the parameter space and that we want the border to be gridlines. (Everything that I say from here on out applies to both x and y, but I'll only state it for x.) For this to be possible, it is necessary and sufficient that N >= n and irem(N-1, n-1) = 0. The condition in the coloring procedure V should be

if irem(round((x-a)*(N-1)/(b-a)), (N-1)/(n-1)) = 0 ...

I haven't tested this extensively, just a little; I just worked it out, literally, on the back of an envelope (actually, two envelopes). So, I wouldn't mind being proven wrong about this. If it matters, I'm using Maple 2016. The reason for all this is that there are only a discrete and evenly spaced set of N values of x that will be passed to V. My formulas are based on knowing exactly what those values are.

Your example plus your statement that the number of visible lines can be varied independently of the grid option make it seem as if the whole continuum of x-values -2..2 are candidates to be passed to the coloring procedure, kind of like how 2-D plot can choose x-values adaptively. Perhaps this is a new feature of Maple 2018's plot3d? If true, it's astounding.

Once again, I think that your idea is superb and very innovative; and I give it a vote up despite its flaws.

@Art Kalb I don't think that parameterizing each gridline would be particularly slow, but there's no point in doing it because you'll run into the same problem that the lines blend into the surface.

@tomleslie You're right, I missed that the "0 to" was missing.

@Art Kalb One solution is as you've already discovered: Use an offset. In particular, you can multiply one dimension by 0.99 in the patchnogrid plot.

Another is to set transparency= 0.1 (or so) for the patchnogrid plot. However, this sometimes creates plots which upon rotation you'll see couldn't possibly exist as objects in three dimensions.

The thickness option applies to the grid lines. Try thickness= 2.