@J4James Sorry, I don't have free access to the article. I am sending you email so that you can send me the PDF attached to an email.

@JohnS Thanks for the compliment. I just added some comments to the above to try to explain exactly what the code is doing.

@J4James Darn it, the link to the PDF doesn't work. Not your fault; it's the "upgrade". Attached files seem to work intermittently. Do you have a URL for the PDF?

@J4James Do you know that the parameters Pr, c, and N are the same for the other plot?

@J4James Okay, I understand what you did to get Nux. And, as far as I can tell, you made no mistakes in transcribing the old system to the new. Now, what is the plot that you are comparing it to? I know that you posted it before, but I can't find it. It may have been a victim of the recent "upgrade". So please post it again.

@awass The multiple-email bug has been fixed.

What you want should be easily achievable by some small modifications to procedure `print/Unit`. Simply removing the square brackets is as easy as

`print/Unit`:= ()-> args;

Adding the color and uprightness should be possible with the Typesetting package, but I don't know the details. I don't know about changing the spacing.

One might as well "endcap" the list right in the procedure:

SignedArea := proc(Pts::list([realcons,realcons]))
local i, P:= [Pts[], Pts[1]];
     add(P[i][1]*P[i+1][2]-P[i+1][1]*P[i][2], i=1..nops(Pts))/2;
end proc:

Now, the bigger challenge is How to detect nonsimplicity of the polygon?

@peardrop How can they be considered linear equations if the a, b, c, d are functions of x[1], ..., x[60]?

@cmcjas Perhaps this will help: i goes from 1 to n. Plugging i = 1 into (i-1)*dxj+a yields a, which is the leftmost point in the first (leftmost) subinterval. Plugging i = n into i*dxj+a yields n*(b-a)/n+a = b -a + a = b, which is the rightmost point in the last (rightmost) subinterval. How's that?

Do you mean that you want to construct a surface from 27 points using c1 = x coordinate, c2 = y coordinate, and c3 = z coordinate?

@J4James No, I don't know how to handle this. In particular, I don't know how you got sys from {Eq1, Eq2, Eq3}. Where does Nux come from? I wonder if there was some mistake made in this process.

What do mean that they are functions of x? Is x distinct from x[1], x[2], ..., x[60]?

@maplelearner

Exactly what plot "shows that the integartion is possible"? My plot (which takes a few minutes to generate) shows that the oscillations do not decrease in amplitude.

plot(R-> Int(r*BesselJ(1,r)*BesselJ(0,r), r= 0..R, digits=4), 100..110);

Also, consider the first term of the asymptotic expansion of the integrand. It does not go to 0 as r -> infinity.

asympt(r*BesselJ(1,r)*BesselJ(0,r), r);

A function is inceasing where its derivative is positive. Typically one finds where the derivative is 0. The zeros divide the real line into intervals. Taking a test point from each interval, determine whether the derivative is positive for that interval.

I can't give more help than that without seeing your function.