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.

@candy898 For positive a, the integral doesn't converge because of the singularity at x=a. You can enter the integral into Maple as

int(ln(x)/(x-a)/(x-1), x= 0..infinity);

Note that it is ln(x), not In(x).

I don't know exactly what Mathematica is doing here. I think it's called Wynn extrapolation or Wynn's epsilon method. Note that Mathematica does express skepticism about its own answers. Also note that the book that you referenced merely explains why Mathematica gives numeric answers for such integrals; it does not claim that those answers are correct (as far as I can tell). I think that the author is more-or-less apologizing for Mathematica's behaviour.

@adel-00 I don't really understand your comment. Are you saying that the program doesn't work for values of theta and phi other than zero? What exactly happens?

@Kitonum Actually, there is no missing parenthesis. The sqrt definitely encompasses the whole expression.

@Kitonum There is a missing right parenthesis in the original. My answer assumed that it closed the whole expression. Yours assumes that it closes the sqrt.

@Adri van der Meer The dynamic Vector approach is significantly faster for sequential output. My earlier Reply was wrong about that.