The advantage is that the intersection points are more clearly separated ==>  better visibility for higher number of lines.

UniformCuts:=proc(N::posint)
local NN:=floor(N/2)*2+1,
      R:=1.25/sin(Pi/2./NN), d:=arccos(1/R),k,p,a,b;
uses plots,plottools;
for k from 0 to N-1 do
  a[k]:=[R*cos(2*k/NN*Pi+d), R*sin(2*k/NN*Pi+d)];
  b[k]:=[R*cos(2*k/NN*Pi-d), R*sin(2*k/NN*Pi-d)];
od:
p:=seq(line(a[k],b[k]),k=0..N-1):
display(disk([0,0],R,color=yellow),p, axes=none, color=red);
end:

UniformCuts(9);

# animation
A:=seq( plots:-display(seq(op(i,%),i=1..k)), k=1..9+1 ):
plots:-display(A, axes=none, insequence, color=red);

Note first that the problem is very ill-conditioned (see previous answer), so unless your coefficients are (almost) exact, the problem is (almost) nonsense.
So, let us assume the coefficients are exact and work with 200 digits!

Unfortunately fsolve refuses to work if the coefficients are complex (but it should!).
Fortunately, using a change of variable it is possible to get real coeffs (in this case) and we can find all the roots via fsolve:

restart;
Digits:=200:
aq := (2.626145*10^(-111)*I)*beta^41+(2.723460372*10^(-55)*I)*beta^25-1.125718*10^(-103)*beta^38-4.42696*10^(-96)*beta^36+(4.038976*10^(-119)*I)*beta^43+(1.897840*10^(-135)*I)*beta^47-1.4537*10^(-128)*beta^44-2.393897*10^(-75)*beta^30-5.345113*10^(-69)*beta^28-(8.88232*10^(-14)*I)*beta^7+(3.78162*10^(-127)*I)*beta^45+(1.87236321*10^(-39)*I)*beta^19-4.22943*10^(-111)*beta^40-5.98764*10^(-82)*beta^32+(1.73215*10^(-27)*I)*beta^13+(3.14576*10^(-144)*I)*beta^49+1.0000002*10^(-150)*beta^50+1.000483*10^(-142)*beta^48+(0.5707492e-4*I)*beta-1.860356732*10^(-56)*beta^26-1.202764308*10^(-50)*beta^24+0.1078870970e-3*beta^4-.1337634356*beta^2+(4.558807*10^(-82)*I)*beta^33+11.99907662+(1.552482870*10^(-49)*I)*beta^23+1.50289073*10^(-33)*beta^16+(4.738072*10^(-89)*I)*beta^35-2.560992731*10^(-45)*beta^22-1.821396189*10^(-40)*beta^20+1.*10^(-159)*beta^52-1.025045*10^(-118)*beta^42-2.708307310*10^(-27)*beta^14-2.880*10^(-137)*beta^46+3.616579272*10^(-22)*beta^12-2.775313578*10^(-17)*beta^10+(2.247453*10^(-75)*I)*beta^31+(4.317773*10^(-69)*I)*beta^29-(8.4742656*10^(-63)*I)*beta^27+(1.086756*10^(-103)*I)*beta^39+(2.875650*10^(-96)*I)*beta^37+(2.927689932*10^(-44)*I)*beta^21+(5.19084*10^(-18)*I)*beta^9+(3.3203077*10^(-35)*I)*beta^17-1.867365177*10^(-8)*beta^6+1.091287414*10^(-12)*beta^8-3.549248092*10^(-36)*beta^18-6.89128*10^(-89)*beta^34-(1.32011*10^(-22)*I)*beta^11+(8.67973*10^(-32)*I)*beta^15+(5.131768*10^(-10)*I)*beta^5-(6.362604*10^(-7)*I)*beta^3-(5.75387*10^(-153)*I)*beta^51:
f:=subs(beta=I*x,aq):
X:=fsolve(f,complex):
B:=X*~I:  # the solutions
seq( evalf(abs(eval(aq, beta=u))), u=[B]): #check
evalf[10]([%]);

[4.183842838*10^(-147), 4.183842838*10^(-147), 6.717604704*10^(-151), 6.717604704*10^(-151), 1.671175541*10^(-152), 1.671175541*10^(-152), 1.531062653*10^(-155), 1.531062653*10^(-155), 2.747698848*10^(-160), 2.747698848*10^(-160), 2.743673043*10^(-162), 2.743673043*10^(-162), 4.127466023*10^(-185), 4.127466023*10^(-185), 8.756351146*10^(-191), 8.756351146*10^(-191), 1.390759103*10^(-192), 1.390759103*10^(-192), 6.530643371*10^(-123), 6.530643371*10^(-123), 1.269161412*10^(-194), 1.269161412*10^(-194), 2.849455036*10^(-195), 2.849455036*10^(-195), 6.407673697*10^(-198), 6.407673697*10^(-198), 7.870992411*10^(-199), 7.870992411*10^(-199), 2.957280721*10^(-198), 2.957280721*10^(-198), 7.563118735*10^(-195), 7.563118735*10^(-195), 3.655881493*10^(-194), 3.655881493*10^(-194), 2.802759752*10^(-192), 2.802759752*10^(-192), 2.265640916*10^(-190), 2.265640916*10^(-190), 1.052966328*10^(-184), 1.052966328*10^(-184), 7.249928057*10^(-162), 7.249928057*10^(-162), 1.960460423*10^(-158), 1.960460423*10^(-158), 1.691321460*10^(-153), 1.691321460*10^(-153), 3.687720329*10^(-151), 3.687720329*10^(-151), 7.453401398*10^(-149), 7.453401398*10^(-149), 1.172077955*10^(-146), 4.492662674*10^(-8)]

The sequence of all the roots is in B.

vx:=[1, 3, 5, 10, 15];
                       [1, 3, 5, 10, 15]
vy:=[2.58, 8.57, 15, 30.10, 45];
                  [2.58, 8.57, 15, 30.10, 45]
f:=Statistics:-LinearFit([1,x], vx, vy, x);
      -0.3817701863353998 + 3.034083850931677 x
p1:=plot(vx,vy, style=point,symbolsize=20,color=red):
p2:=plot(f,x=0..16):
plots:-display(p1,p2,gridlines);

So, you want to generate an n x n random unitary matrix, k (<n) columns being fixed.
You could do the following:
1. Generate an n x k matrix A with the fixed columns.
2. Generate an  n x (n-k) random matrix B.
3. C := <A | B>   (=square matrix)
4. Check that det(C)<>0; if not, goto 2.
5. Apply the Gram-Schmidt process to C resulting a matrix G; it will be unitary.
6. Move the first k columns of G to their original positions. This is the desired matrix.
 

with(plots):  with(plottools):
 H:=(20/3)*sqrt(3.):
 Sph := sphere([0, 0, 0], 10,   color=gold, transparency=0.4):
 Cyls := display([seq(cylinder([0,0,h], sqrt(100-h^2),2*abs(h), color=red, strips=50),
     h=-10 .. -H/2.0, 0.2)], insequence):
 display([Sph, Cyls], style=surface, lightmodel=light1, axes=none);

dsolve({deqv, v(0) = v__0}, v(s)) assuming k>0, v__0>0;

A space filling curve very easy to implement (and understand) is due to Lebesgue and simplified by Schoenberg:

F:=proc(t)
   local u:=abs(2*round(t/2)-t);
   piecewise(u<=1/3,0, u<2/3, 3*u-1,1)
end:
PN:=4:
PX:=proc(t)  add(F(9^k*t)/2^(k+1),  k=0..PN) end:
PY:=proc(t)  add(F(9^k*3*t)/2^(k+1),k=0..PN) end:

plot(PX,0..1,numpoints=1000);

plot([PX,PY,0..1], numpoints=10000);

(Of course the real space filling curve is for PN=infinity).

If you want the numeric result, use:

int(f,-1..1,numeric);
   1.177418042
For a symbolici result, Markiyan suggestion should work but it gives a wrong answer (Maple's bug!):

g := proc(x) local r; r := sin(x)+x*cos(x); if is(abs(r) < 1/2) then sin(x) else cos(x) end if; end proc;

int(g, -1..1);evalf(%);
    2*sin(1)
   1.682941970

A piecewise attempt also fails for symbolics!

Workaround for a symbolic solution:

x0:=RootOf(2*_Z*cos(_Z)+2*sin(_Z)-1,0..1):
int(cos(x),x=-1..-x0)+int(sin(x),x=-x0..x0)+int(cos(x),x=x0..1);
     2*sin(1)-2*sin(RootOf(2*_Z*cos(_Z)+2*sin(_Z)-1, 0 .. 1))
evalf(%);
    1.177416231

EDIT. The value of the minimum given by minimize is of course wrong (bug!).

loc:=%[2][][1];
               {x = 3405.798701, y = 2169.889609}
eval([eq1,eq2,eq3],loc);
                       [0.01, 0., -0.01]

P.S. You may be interested in  ?geometry,RadicalCenter

EQ:=convert(lhs(eq1),set) union convert(lhs(eq2),set):
sol:=solve(EQ):
eval([fn1,gn1],sol);

N.B. You should avoid evalm and use "+", ".".

Use:

[A$2] .~ LM;

or

(A$2) .~ LM;

GetCoeff:=proc(f,u)
  local x:=indets(f,name);
  if u=0 or u=1 then return eval(f,x=~0) fi; 
  eval(foldl(coeff,f,op(u)), x=~0);
end proc:

It is possible to represent (and "see")  only surfaces. You cannot see in the interior of a 3-dimensional set in R^3 (between any interior point and your eye there exist an infinity of other interior points obfuscating the first one).

Today you have problems easier by hand than with Maple :-)

f = (x+y)/(x+y+z), g = (x+z)/(x+y+z), h = (y+z)/(x+y+z)  have the same integral (by symmetry).

==> int(f) = 1/3*int(f+g+h) = 1/3*int(2) = 1/3*2*(4*Pi/3)/8 = Pi/9.

Edit.
A real challenge for any CAS would be to compute the integral of
(x^2017+y^2017)/(x^2017+y^2017+z^2017)
in the first octant of the unit ball. [=Pi/9]

dxdy  :=   y=-sqrt(1-x^2) .. sqrt(1-x^2), x=-1..1 ;

int(sin(x^2) *cos(x^2) , dxdy) =

1/2*int( sin(x^2+y^2) + sin(x^2 - y^2), dxdy ) =

1/2*int( sin(x^2+y^2) , dxdy)  + 0 =

1/2* int( sin(r^2)*r, r=0..1, t=0..2*Pi ) =

1/2*(1 - cos(1))*Pi =

Pi * sin(1/2)^2