@vv I didn't mean to imply that extending to rational functions would involve simply changing A to rational functions. You need to multiply each row by the LCM (or just the product) of its denominators, make an estimate for the upper bound of the degree of the determinant of the new matrix (you can use the sum of the max of the degrees in each row), compute the determinant by the present process (with the number of imaging points determined by that degree estimate), and divide the result by the product of the LCMs.

@Rouben Rostamian  Yet the following shows how to using imaging and interpolation to compute the determinant of a 170 x 170 dense matrix of linear univariate polynomials with integer coefficients in 40 seconds:

restart:
n:= 170:
A:= LinearAlgebra:-RandomMatrix((n$2), generator= (()-> randpoly(x, degree= 1))):
Af:= (x::integer)-> rtable(eval(A, :-x= x), datatype= integer):
DET1:= CodeTools:-Usage(
   CurveFitting:-PolynomialInterpolation(
      [seq([k,LinearAlgebra:-Determinant(Af(k))], k= -iquo(n,2)..iquo(n,2)+1)],
      x
   )
):

memory used=0.82GiB, alloc change=48.03MiB, cpu time=39.88s, real time=40.29s, gc time=218.75ms

This can be easily distributed over multiple cores, like this (my present machine has 4 cores):

save DET1, `DET1.m`: #So we can compare results.
restart:
n:= 170:
A:= LinearAlgebra:-RandomMatrix((n$2), generator= (()-> randpoly(x, degree= 1))):
DET2:= CodeTools:-Usage(
   CurveFitting:-PolynomialInterpolation(
      [Grid:-Seq['tasksize'= 1](
         [k, LinearAlgebra:-Determinant(rtable(eval(A, x= k), datatype= integer))], 
         k= -iquo(n,2)..iquo(n,2)+1)
      ], x
   )
):

memory used=23.61MiB, alloc change=0 bytes, cpu time=1.02s, real time=18.02s, gc time=625.00ms

(That "memory used" and "cpu time" must only be measured for the master process, but the "real time" is real.)

read `DET1.m`:
evalb(DET1 = DET2);

     true

Extending this idea to A's entries being univariate rational functions with rational coefficients is not significantly more complicated.

This is essentially Rouben's method, but using only one line of code after the solve:

((F,B)-> eval(B, lhs~(F) =~ 0))(selectremove(evalb, sol_full));

Please post that code in plaintext so that I can copy-and-paste it. If you copy-and-paste directly from Maple to MaplePrimes, and you don't select "as image", that should do it.

@grzybs It's unfortunate that the conversion of a definite iterated integral to a different coordinate system requires several meticulous steps.

@Ronan Because of the different ways that mixtures of Pi and floats are treated in different versions of Maple, the Pi should be inside the brackets of the evalf. And you might as well put the 180 in there also. So,

evalf(FindAngle(p1,p2)*180/Pi);

@grzybs It looks like you forgot to copy-and-paste your new values, and simply retained Kitonum's values in your Reply above. Anyway, have a look at my worksheet below. The only thing that you need to change is the second-to-last command, which gives values to the constants.

So, did you find the above model useful or not? Some feedback would be appreciated.

@BlueSourBalls Your plot is totally correct, and your code is very good for a beginner. One tip: Don't use square brackets in place of parentheses. You have a stray square bracket in the line beginning f(x):=.

Since you indicate that you have a graph, please show the commands that you used to obtain it. Or show the graph.

@taro I don't see how it could possibly be even close to as good as MaplePrimes. There's already a huge unused question-answering capacity here.

@taro The plane is tangent to the sphere. The intersection is the single point (1,1,1); the intersection is not a circle.

What makes you think that your ability to divide the region into subintervals produces results any better (faster, more accuate, etc.) than Maple's own ability to determine how to do the integration over the whole domain?

You use epsilon = 1e-3, so the error in each subinterval can be as much as 1 part in 1000. Then you add together 5000+ such results. So, the end result has 0 significant digits (5000 x 1/1000 > 1).

@dharr Your technique works except for when $ is used as the end-of-parameters marker. In that case, ` $` will appear in op(1, eval(f)). Note that that's ` $`, not `$` (there's a space at the beginning).

Also, there's a distinction between parameters  and arguments, which I elaborate upon in my Answer below.

@ Yes, I explicitly wrote it to retain any trailing zeros that were included when the number was entered. I can remove them if you want. Trailing zeros convey information about precision; leading zeros (to the left of the decimal point) do not.