Here is an independent solution of the logical problem. Code is quite complicated, because based on the phase-solving, but it works fast. Carl, do  you get the same answers in your package?

restart;

Position:=[$ 1..5]:

Names:= [Bob, Keeley, Rachael, Eilish, Amy]:

TV:=[Simpsons, Coronation, Eastenders, Desperate, Neighbours]:

Dest:= [Fra, Aus, Eng, Afr,Ita]:

Ages:= [14, 21, 46, 52, 81]:

Hair:=[afro, long, straight, curly , bald]:

Lives:= [town, city, village, farm, youth]:

P:=proc(A,B)

if (A implies B) and (B implies A) then true else false fi;

end:

L:=[]:

for i in Position do

for j in Names do

for k in TV do

for p in Dest do

for q in Ages do

for s in Hair do

for t in Lives do

if convert([P(i=3,k=Desperate), P(i=1,j=Bob), (k=Simpsons implies t<>youth) and (t=youth implies k<>Simpsons), (p=Afr implies j<>Rachael) and (j=Rachael implies p<>Afr), P(t=village,q=52), P(p=Aus,s=straight), P(p=Afr,k=Desperate), P(q=14,i=5), P(j=Amy, k=Eastenders), P(p=Ita,s=long), P(j=Keeley,t=village), P(q=46,s=bald), P(i=4, p=Eng),  (k=Coronation implies s<>afro) and (s=afro implies k<>Coronation), (j=Rachael implies s<>afro) and (s=afro implies j<>Rachael), P(q=21,t=youth), P(k=Coronation, s=long), P(q=81,t=farm), P(p=Fra,t=town), (j=Eilish implies s<>straight) and (s=straight implies j<>Eilish)], `and`) then L:=[op(L), [i,j,k,p,q,s,t]] fi:

od: od: od: od: od: od: od:

M:=[ListTools[Categorize]((x,y)->x[1]=y[1], L)]:

R:=[]:

for i in M[1] do

for j in M[2] do

for k in M[3] do

for p in M[4] do

for q in M[5] do

if ((i[3]=Simpsons and j[7]=youth) or (j[3]=Simpsons and (i[7]=youth or k[7]=youth)) or (k[3]=Simpsons and (j[7]=youth or p[7]=youth)) or (p[3]=Simpsons and (k[7]=youth or q[7]=youth)) or (q[3]=Simpsons and p[7]=youth)) and ((i[2]=Rachael and (j[4]=Afr or k[4]=Afr or p[4]=Afr or q[4]=Afr)) or (j[2]=Rachael and (k[4]=Afr or p[4]=Afr or q[4]=Afr)) or (k[2]=Rachael and (p[4]=Afr or q[4]=Afr)) or (p[2]=Rachael and q[4]=Afr)) and ((i[3]=Desperate and j[3]=Neighbours) or (j[3]=Desperate and (i[3]=Neighbours or k[3]=Neighbours)) or (k[3]=Desperate and (j[3]=Neighbours or p[3]=Neighbours)) or (p[3]=Desperate and (k[3]=Neighbours or q[3]=Neighbours)) or (q[3]=Desperate and p[3]=Neighbours)) and ((i[3]=Coronation and j[6]=afro) or (j[3]=Coronation and (i[6]=afro or k[6]=afro)) or (k[3]=Coronation and (j[6]=afro or p[6]=afro)) or (p[3]=Coronation and (k[6]=afro or q[6]=afro)) or (q[3]=Coronation and p[6]=afro)) and ((i[2]=Rachael and j[6]=afro) or (j[2]=Rachael and (i[6]=afro or k[6]=afro)) or (k[2]=Rachael and (j[6]=afro or p[6]=afro)) or (p[2]=Rachael and (k[6]=afro or q[6]=afro)) or (q[2]=Rachael and p[6]=afro)) and not ((i[2]=Eilish and j[6]=straight) or (j[2]=Eilish and (i[6]=straight or k[6]=straight)) or (k[2]=Eilish and (j[6]=straight or p[6]=straight)) or (p[2]=Eilish and (k[6]=straight or q[6]=straight)) or (q[2]=Eilish and p[6]=straight))  then

R:=[op(R), [i,j,k,p,q]] fi:

od: od: od: od: od:

Result:=select(x->convert([seq(nops({seq(x[i,j], i=1..5)})=5, j=1..7)],`and`), R):

print(`First solution`);

Matrix(Result[1])^+;

print(``);

print(`Second solution`);

Matrix(Result[2])^+; 

Puzzle.mws

with(LinearAlgebra):

a, b, c:=<22,11,5>, <13,6,3>, <-5,-2,-1>:

a1:=a/norm(a,2):

b0:=b-(a1.b)*a1:

b1:=b0/norm(b0,2):

c0:=c-((a1.c)*a1+(b1.c)*b1):

c1:=c0/norm(c0,2):

a1,  b1,  c1;

Vectors  a1, b1, c1  are unit and mutually perpendicular:

a1.a1, b1.b1, c1.c1, a1.b1, a1.c1, b1.c1;

            1, 1, 1, 0, 0, 0

convert(op(0.25)[1], symbol)/(10^(-op(0.25)[2]));

 ``(op(0.25)[1])/``(10^(-op(0.25)[2]));

The second method allows you to use this fraction in subsequent calculations:

expand(%*3/5);

              3/20

Solution without  dsolve (only for rectangular tables):

restart;

A:=plots[animate](plot,[[2/Pi*abs(arcsin(sin(Pi/2*(0.5+0.2*t)))), 2/Pi*abs(arcsin(sin(Pi/2*(0.5+sqrt(3)/5*t)))), t=0..a], thickness=5, color=red, numpoints=1000], a=0..20, frames=200):

B:=plottools[rectangle]([0,1], [1,0], thickness=2, color=green):

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

with(RealDomain):

 plot(x^(1/3), x=-8 .. 8, thickness=5, scaling=constrained);

For Classic M 12:

A:=1/2+(1/3)*sqrt(2);

numer(A)/(``*denom(A));

For Standard:

A := 1/2+(1/3)*sqrt(2);

numer(A)/convert(denom(A), symbol);

Another variant:

A := 1/2+(1/3)*sqrt(2);

``(numer(A))/denom(A);

For your example it's more reliable to use a linear model, which is obtained if to apply  ln  function to your equation:

ln(x)=ln(a)+b*ln(y)+c*ln(z)  or  x1=a1+b*y1+c*z1 , where  x1=ln(x) ,  a1=ln(a) ,  y1=ln(y) ,  z1=ln(z) 

xlist := [0.1e-1, 0.1e-1, 0.1e-1, 0.113e-2, 0.144e-2, 0.171e-2, 0.21e-2, 0.23e-2, 0.183e-2, 0.233e-2, 0.28e-2, 0.331e-2, 0.409e-2, 0.304e-2, 0.377e-2, 0.459e-2, 0.576e-2, 0.708e-2, 0.424e-2, 0.502e-2, 0.623e-2, 0.763e-2, 0.918e-2, 0.529e-2, 0.63e-2, 0.778e-2, 0.965e-2, 0.116e-1, 0.63e-2, 0.747e-2, 0.899e-2, 0.113e-1, 0.137e-1, 0.724e-2, 0.872e-2, 0.105e-1, 0.13e-1, 0.158e-1, 0.848e-2, 0.1e-1, 0.122e-1, 0.15e-1, 0.182e-1]:

ylist := [16.529, 13.333, 10.929, 2.66, 2.75, 2.77, 2.69, 2.58, 4.4, 4.67, 4.62, 4.51, 4.48, 7.5, 7.58, 7.64, 7.75, 7.77, 10.1, 10.1, 10.1, 10.2, 10.2, 12.6, 12.7, 12.8, 12.7, 12.8, 14.9, 15, 15, 15, 15, 17.3, 17.5, 17.4, 17.3, 17.3, 20, 19.9, 19.9, 19.9, 20]:

zlist := [100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000, 1, 10, 100, 1000, 10000]: 

X:=convert(map(ln,xlist), Vector):

Y:=convert(map(ln,ylist), Vector):

Z:=convert(map(ln,zlist), Vector):

YZ:=Matrix([Y,Z]):

Statistics[LinearFit]([1,y1,z1], YZ, X, [y1,z1]):

a=exp(tcoeff(%)), b=coeff(%,y1), c=coeff(%,z1);

            a = 0.0004189912543,  b =0 .996002908082386118, c = 0.0849219222885009006

L1:=x + y - 1 = 0, x - z - 1 = 0:

L2:=x + y - 7 = 0, x - y + 1 = 0:

A:=plots[intersectplot](L1, x=-5..5, y=-5..5, z=-3..5, thickness=2, color=red):

B:=plots[intersectplot](L2, x=-5..5, y=-5..5, z=-3..5, thickness=2, color=blue):

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

For your data, you can find exactly unknown a, b, c:

X:=[0.01, 0.01, 0.01]:

Y:=[16.526, 13.333, 10.929]:

Z:=[100, 1000, 10000]:

f:=x=a*y^b*z^c:

solve({seq(eval(f, {x=X[i], y=Y[i], z=Z[i]}), i=1..3)});

                             {a = 0.01, b = 0., c = 0.}

Let the numbers assigned to each chair from 1 to 8. 

Seating for 8 people means to establish a one-to-one correspondence between these people and these numbers. Total we get  96 solutions, including 6 unique solutions:

restart;

Names := [Alan, Amy, Brad, Beth, Charles, Debbie, Emily, Frances]:

P := proc (n1, n2)

`if`(irem(abs(n2-n1), 7) <= 1, true, false) end proc:

Permutations := combinat[permute]([`$`(1 .. 8)]):

Arrangements := []:

for p in Permutations do

if eval(`and`(`and`(`and`(`and`(`and`(`and`(`and`(`and`(P(Amy, Alan), P(Brad, Beth)), `or`(P(Charles, Debbie), P(Charles, Emily))), P(Frances, Debbie)), `and`(`and`(`and`(`not`(P(Amy, Brad)), `not`(P(Amy, Beth))), `not`(P(Alan, Brad))), `not`(P(Alan, Beth)))), `and`(`not`(P(Brad, Charles)), `not`(P(Brad, Frances)))), `not`(P(Debbie, Emily))), `and`(`not`(P(Alan, Debbie)), `not`(P(Alan, Emily)))), `not`(P(Amy, Charles))), {seq(Names[i] = p[i], i = 1 .. 8)})then Arrangements := [op(Arrangements), p] end if end do:

seq([seq([Names[i], p[i]], i = 1 .. 8)], p = Arrangements);

Arrangement_Birthda.mw

restart;

plot(x^2, x=0..2, thickness=2, view=[-2..2, -0.5..4.5], scaling=constrained);

Lorny! Your system can be solved exactly with respect to  T  and  C  at any value of  n2 . To do this,  n1  and  n3  should be defined exactly as a fractions. Solutions are expressed in terms of the roots of cubic equations and even can be expressed by radicals, but very cumbersome expressions. Simple expression have only solutions for  n2 = 3/10 . If  n2 <0.3  or  n2> 0.3246543748 (approximately), the system has a unique solution, if  n2 = 0.3  or  n2 = 0.3246543748 (approximately), then two solutions, and if  0.3 <n2 <0.3246543748, then three solutions:

n1:=1: n3:=3/10:

n2:=1/5:

solve({T= n1/(1+abs(T-S)), S=n2/(n3+abs(T-S))});

n2:=31/100:

solve({T= n1/(1+abs(T-S)), S=n2/(n3+abs(T-S))});

Several specific solutions:

n1:=1: n3:=0.3:

n2:=0.29:

solve({T= n1/(1+abs(T-S)), S=n2/(n3+abs(T-S))});

n2:=3/10:

solve({T= n1/(1+abs(T-S)), S=n2/(3/10+abs(T-S))});

n2:=0.31:

solve({T= n1/(1+abs(T-S)), S=n2/(n3+abs(T-S))});

n2:=0.3246543748:

solve({T= n1/(1+abs(T-S)), S=n2/(n3+abs(T-S))});

n2:=0.3246543749:

solve({T= n1/(1+abs(T-S)), S=n2/(n3+abs(T-S))});

The plot of solutions. The plot, where the solution is not unique, isolated vertical dashed lines:

> restart;

n1:=1: n3:=0.3:

f:=n2->solve({T= n1/(1+abs(T-S)), S=n2/(n3+abs(T-S))}):

A:=plot(['eval(T, f(n2))', 'eval(S, f(n2))'], n2=0..0.2999, color=[red, blue], thickness=2):

B:=plot(['eval(T, f(n2))', 'eval(S, f(n2))'], n2=0.3248..1, color=[red, blue], thickness=2):

E:=plot(['eval(T, f(n2)[1])', 'eval(T, f(n2)[2])', 'eval(T, f(n2)[3])'], n2=0.3001..0.3246, color=red, thickness=2):

F:=plot(['eval(S, f(n2)[1])', 'eval(S, f(n2)[2])', 'eval(S, f(n2)[3])'], n2=0.3001..0.3246, color=blue, thickness=2):

G:=plot([0.3, t, t=0..1], linestyle=3, color=black):

H:=plot([0.3247, t, t=0..1.001], linestyle=3, color=black):

plots[display](A, B, E, F, G, H, view=[0.25..0.4, 0.4..1.2]);

Easiest to use  add  function:

PartialSum:=proc(n::posint)

evalf(add(4*(-1)^(k-1)/(2*k-1), k=1..n));

end proc:

Example:

PartialSum(1000);

            3.140592654

To solve the second problem, we use property an alternating series: the remainder of a series does not exceed the first of the neglected terms (in absolute value).

solve(4/(2*n-1)<10^(-5));

Therefore, we take n = 200001

PartialSum(200001);

                      3.141597654

evalf(Pi - %);

                      -0.5*10^(-5)

n1:=1: n3:=0.3:

f:=n2->fsolve({T= n1/(1+abs(T-S)), S=n2/(n3+abs(T-S))}):

plot(['eval(T, f(n2))', 'eval(S, f(n2))'], n2=0.5..3, color=[red, blue], numpoints=1000);

Replace  e^(-1.5*t)  by  exp(-1.5*t) :

A:=Int(exp(-1.5*t)/sqrt(t*(t+1)), t = 1 .. infinity);

evalf(A);