@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.

@valery12 Your code for 4th-order Runge-Kutta is already correct! You just need to specify the differential equation. That is, you need to set the value of f!

To clarify (possibly for other readers), you want to increase the plot's extent in the second coordinate dimension, which is not quite the same thing as the y-coordinate because you used y as a parameter.

In your Question, you say z/(1+z), but in the attached code, you use z/(1-z). Which is correct?

@valery12 What Mac Dude means is that you didn't set the values of f or _Meth.

@briceM No, you can't use cmaple to run .mw files (Maple worksheets). You need to extract the code from the worksheets into plaintext files first.

That's certainly not C code. I think that it's Fortran-ish. It's not quite Fortran because the comments are not distinguished from the code in any obvious way. Before we try translating this to Maple, could you supply code that's known to work in its source language?

What you describe is possible, but there's likely to be a much better way: saving the entries or the whole Jacobian as procedure code, possibly compiled, which may take much less memory because of the encapsulation of common subexpressions. That encapsulation can be done automatically.

Before we proceed, I need to know what you need to do with the Jacobian afterwards. Will you just use it for numeric evaluation? or will you use it for symbolic analysis?

@Adam Ledger SCR = Software Change Request