You specify an undirected graph, so the distance matrix should be symmetric. But in your graph, for example,  (2, 22) = 547  and  (22, 2) = 500

You can easily find the general solution of this equation:

restart;

g :=0.88641:

e := 2.53128:

eq := tan(g)-e*sin(f)/(1+e*cos(f)):

solve(eq, AllSolutions);

                    -2.566269006+6.283185307*_Z1,   1.197496352+6.283185307*_Z1

about(_Z1);

                                       Originally _Z1, renamed _Z1~:
                                          is assumed to be: integer

restart;

a := 0:  b := 2:

n[x] := 20:  n[t] := 50:  dt := 0.01:  c := 1:  dx := (b-a)/n[x]:

for i while i <= n[x] do

  if 5 <= i and i <= 10 then u[i] := 2 else u[i] := 1:

  end if;

 end do;

seq(u[i], i = 1 .. 20);

                                    1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1

P:=proc(N::posint, F::procedure)

uses plottools;

seq(seq(disk([i,j], 0.1, color=`if`(F(i,j), red, blue)), i=1..N), j=1..N);

end:

Example:

plots[display](P(7, (i, j)->`if`(i>j, true, false)), axes=none);

 Improvement.

 Instead

plots[display](P(7, (i, j)->`if`(i>j, true, false)), axes=none); 

it is shorter to write 

plots[display](P(7, (i, j)->is(i>j)), axes=none);

plots[implicitplot]([seq(x^2-2*y^2=C, C=-5..5)], x= -5..5, y= -5..5, numpoints=10000, scaling=constrained);

Numpoints  option allows you to significantly improve the quality of construction, and  scaling=constrained  option equalizes the scale of the axes.

Eq1:=diff(f(x), x)=1, g(x)=diff(f(x), x);

cs:=f(0)=0;

p:=dsolve({Eq1, cs}, numeric);

plots[odeplot](p, [x, g(x)], 0..10);

If you open the brackets in the canonical equation of a circle  (x-x0)^2+(y-y0)^2=R^2 , you get an equation of the form  x^2+y^2+a*x+b*y+c=0 . 

Your problem can be easily solved symbolically, ie exactly. Here is an example:

restart;

F := x^2+y^2+a*x+b*y+c:

XY := [[0, 0], [2, 3/2], [3, 1], [2, -1]]:

minimize(add((eval(F, {x = XY[i, 1], y = XY[i, 2]}))^2, i = 1 .. nops(XY)), location);

assign(op(%[2])[1]):

A := plot(XY, style = point, symbolsize = 15, symbol = solidcircle, color = red):

B := plots[implicitplot](F, x = -1 .. 4, y = -2 .. 2, color = green):

plots[display](A, B, scaling = constrained);

Parameters  x0, y0, R  can be easily expressed in terms of  a, b, c . This you can make yourself.

plots[implicitplot](0.3939708949*x^2 - 0.005975799853 + 0.6345432059*y^2=0,

x=-1..1, y=-1..1, thickness=2, numpoints=50000, scaling=constrained);

Addition.  We can do the job without  numpoints option if to shorten the axis, because the ellipse is very small:

plots[implicitplot](0.3939708949*x^2 - 0.005975799853 + 0.6345432059*y^2=0,

x=-0.2..0.2, y=-0.2..0.2, thickness=2, scaling=constrained);

Irem:=proc(a::integer, b::integer)

b*frac(a/b);

end;

Example:

Irem(8, 3); 

     2

combine(-ln(x)+ln(y))  assuming x>0,y>0;

with(RealDomain):

Sol:=solve({x^2 + 2*y*z^2 = 0, y^2  - 3*x*z = 0, 1/3*x*y^2 + 2*y*z^3 = 0, -1/3*y^4 - 6*y*z^4 = 0}):

L:=allvalues({Sol}):

op({seq(simplify(L[i], power), i = 1 .. 3)});  # Final answer

We got infinitely many solutions, depending on one real parameter  z .

Addition: Geometrically, the set of solutions is the union of two space lines in R^3. First line coincides with the axis Oz.

A:=plots[spacecurve]([0,0,z], thickness=3, color=red, z=-5..5, numpoints=1000):

B:=plots[spacecurve]([1/3*z*(-18*z)^(2/3),(-18*z)^(1/3)*z,z], thickness=3, color=red, z=-5..5, numpoints=1000):

plots[display](A, B, axes=normal, labels=[x,y,z], view=[-5..5, -5..5, -4..4], scaling=constrained); 

Divergence of the integral is easily obtained by a direct calculation using the result by Alejandro:

MultiSeries:-series(sin(x)^cos(x), x = Pi, 2);

limit(int(%, x = 0 .. Pi - t), t = 0, right);

I understood the original question in the sense that it is required to find the sum of all possible values  f(x_1, x_2, .. , x_n) ​​, where each argument  x_i  independently of the others runs from  0  to  k_i .

Summa:=proc(f::procedure, L::list)

local n, T, S;

uses combinat;

n:=nops(L);

T:=cartprod([seq([seq(i, i=0..L[k])], k=1..n)]);

S:=0;

while not T[finished] do S:=S+f(op(T[nextvalue]()))

od;

S;

end; 

Example:

Summa((x1,x2,x3)->f(x1,x2,x3), [2, 3, 5]);

Another example (random polynomial in two variables, where the degree of each variable does not exceed 3, and the coefficients are random integers in the range -9 .. 9):

roll := rand(-9 .. 9):

Summa((a,b)->roll()*x^a*y^b, [3,3]);

The problem is easily solved by hand.

The answer in the interval  0..Pi:   RealRange(0, Open(Pi/6)), Pi/3, 2*Pi/3, Pi

This solution is periodically extended to the whole real axis with period  Pi .

Illustration to the answer:

plot([abs(tan(x)*tan(2*x)*tan(3*x))+abs(tan(x)+tan(2*x))-tan(3*x), [(1/6)*Pi, t, t = -1 .. 20], [(1/4)*Pi, t, t = -1 .. 20], [(1/2)*Pi, t, t = -1 .. 20], [2*Pi*(1/3), t, t = -1 .. 20], [3*Pi*(1/4), t, t = -1 .. 20], [5*Pi*(1/6), t, t = -1 .. 20]], x = 0 .. Pi, -1 .. 20, thickness = [3, 1, 1, 1, 1, 1, 1], linestyle = [1, 3, 3, 3, 3, 3, 3], color = [red, seq(black, i = 1 .. 6)], discont = true, tickmarks = [spacing((1/6)*Pi), default], numpoints = 3000);

Edited: the plot is improved. Dashed lines correspond to the points of discontinuity.

This software package  

http://www.umsolver.com/usa/mathematics/free-algebra-equation-solver/

makes it possible to solve step by step much more problems.