You can simply make a screenshot of a web page  and paste it into a working document Maple or initially to process in Paint (eg cut) and then paste.

In help  plot3d Color Functions  you can read:  "In general, specifying the color using procedures or expressions works only with surfaces defined explicitly by a procedure or expression.  For example, spacecurves generated by the plots[spacecurve] command and most surfaces generated by the plottools package commands cannot be colored in this way.  Instead, you can color those surfaces as described in plot/color."

See this  article  from wiki.

Here are plots of some such 3-ellipses:

P:=[x,y]: A:=[-1,0]: B:=[2,0]: C:=[0,1]:
F:=sqrt(add(`^`~(P-A,2)))+sqrt(add(`^`~(P-B,2)))+sqrt(add(`^`~(P-C,2)));
plots[implicitplot]([F=4, F=5, F=6], x=-3..3, y=-3..3, scaling=constrained, gridrefine=3, color=[red,blue,green], view=[-2..3,-2..2.5]);

In these examples, we see that some foci (points A, B, C) are located outside of the curve.

If you mean the polar equation for any circle, then your equation  r^2=(r*cos(t)-a)^2+(r*sin(t)-b)^2  is wrong, because at the right  r  is the polar radius for a point on the circle and at the right there is the radius of the circle. Here is the correct calculation with a plotting:

solve((r*cos(t)-a)^2+(r*sin(t)-b)^2=R^2, r);
plot(eval([%], {a=-1,b=0, R=2}), t=0..2*Pi, color=red, coords=polar);

The obvious solution for the general matrix - write it down as a procedure. As the input of this procedure can be submitted any specific  nxn matrix:

ModalMatrix:=proc(A)
uses LinearAlgebra;
simplify(fnormal~(evalf~(Eigenvectors(A)[2])), zero);
end proc:

Examples of use:

ModalMatrix(<-4,-2; 3,1>);
ModalMatrix(<-4,-2,5; -2,1,7; 5,7,9>);
                     

Since starting with n = 3 symbolic solution can be very cumbersome, simplifying the commands used.

Simple  Rand  procedure solves your problem. It returns a random number in a range  -r .. r . r  can be specified in any format: integer, fraction, float, symbolic.

Rand:=proc(r::positive)
local rf, p;
rf:=convert(r, float);
p:=rand(-rf .. rf);
p();
end proc:

Example of use ( 10 random numbers in the range  -Pi .. Pi ):

seq(Rand(Pi), i=1..10);

0.065963526, 0.743846604, 0.980407970, 1.155281063, 0.793998054, -3.075735010, 0.546438190, -1.067155770, 2.994200916, 1.581718683

The procedure was written in Maple 2016. In older versions you may need  the call of  RandomTools[Generate]  command.

A bit shorter:

A := simplify([a=log[12](27), log[36](24)]);
B := simplify(subs(ln(3)=k*ln(2), A));
simplify(eval(B[2], k=solve(B[1], k))); 

plot3d([r*cos(theta), r*sin(theta), 1-r^2*cos(theta)^2], r = 0 .. 1, theta = 0 .. 2*Pi, color = khaki, filled = true, style = surface, numpoints = 40000);

The formal parameters of the procedure  P  are the inital values of the velocity  v1  and the coordinate of the point  x1  and the number of the steps  n :

P:=proc(v1, x1, n)
local F, v, x, m, Δt, i;
global V, X;
F:=rand(-0.5..0.5);
v[1]:=v1; x[1]:=x1; m:=1; Δt:=1;
for i from 2 to n do
v[i] := v[i - 1] +(1/m)*F()*Δt;
x[i] := x[i - 1] + 0.5*(v[i - 1] + v[i]);
od;
V:=convert(v,list);
X:=convert(x,list);
V, X;
end proc: 
 

Example of use:

P(1, 0, 10);

[1, .9794824628, .7573440398, 1.053436739, 1.478762216, .9900980562, .4951472162, .6204418912, 1.056397926, .6999520291],   [0, .9897412315, 1.858154483, 2.763544872, 4.029644350, 5.264074486, 6.006697122, 6.564491676, 7.402911584, 8.281086562]

Use  solve  command instead of  fsolve  for your example. Here a simplier example of the matrix  <x,1; 1,x> , which leads to the solution of the system:

solve({-c1+c2=0, c1-c2=0, c1^2+c2^2=1}, explicit);

     {c1 = (1/2)*sqrt(2), c2 = (1/2)*sqrt(2)}, {c1 = -(1/2)*sqrt(2), c2 = -(1/2)*sqrt(2)}

a:=convert(2.1345, string);
b:=[seq(a[i], i=1..length(a))];
remove(t->t=".", b);
parse~(%);

Another way (shorter):

ListTools:-Reverse(convert(op(1, 2.1345), base, 10));
                                              

Error reason is that the syntax of  VectorCalculus:-int  is different from the  int  command  from  Maple kernel, which you have used. Therefore, do not download the whole package  VectorCalculus, and use the long form of command call as  VectorCalculus:-command

I corrected only the first error in the calculation of the volume. The remaining errors you can easily fix yourself.

See the corrected file.

error_integral1.mw

Slightly reduce the right end of  x - range:

f:=x->piecewise(x>-Pi and x<=Pi, 3*x):
plot(f(x), x=-Pi..Pi, discont=true):
a:=-Pi:  b:=Pi:
p:=b-a:
fp:=x->f(x-floor((x-a)/p)*p):
plot(fp(x), x=-6*Pi..7*Pi-0.001, discont=true, tickmarks=[piticks,piticks]);

Example:

A:=plot3d([cos(phi)*sin(theta), sin(phi)*sin(theta), cos(theta)], phi=-arccos(0.9)-0.01 .. Pi/2-arccos(0.9), theta=arccos(0.9) .. `if`(phi>=-arccos(0.9) and phi<=arccos(0.9), arcsin(0.9/cos(phi)), Pi/2), style=surface, color=khaki, scaling=constrained, numpoints=25000):
B:=plots[display](seq(plottools[rotate](plots[display](A), Pi/2*k, [[0,0,0], [0,0,1]]), k=0..3)):
C:=plottools[reflect](B, [[0,0,0], [1,0,0], [0,1,0]]):
Cube:=plottools[cuboid]([-0.9,-0.9,-0.9], [0.9,0.9,0.9], color=pink, transparency=0.8):
plots[display](B, C, Cube, axes=none);

Edit.
 

For example,  if  i=4  then

`or`(seq(R[i]=R[k], k=1..i-1))   means

R[4]=R[1]  or  R[4]=R[2]  or  R[4]=R[3]