Let's consider the vertices of a regular hexagon with a side of 1, the center of which lies at the origin. For definiteness, we will assume that the point in the center is red. Since there are still 6 points left, there are a total of 2^6=64 coloring options (see the figure below). Of these 64 options, only 3 options do not have an equilateral triangle with equally colored vertices. Let's take one of these options (the next figure) and add 3 more points to the right, forming 2 equilateral triangles with a side of 1 with the right red point. If the rightmost point is black, then we immediately get an equilateral triangle with black vertices with a side  sqrt(3). If this point is red, then with any choice of colors for the remaining 2 points we also get an equilateral triangle with vertices of the same color.

restart;
with(plots):
L:=combinat:-permute([red$6, black$6], 6):
P:=seq(pointplot([[0,0],seq([cos(2*Pi*k/6),sin(2*Pi*k/6)], k=0..5)], color=[red,p[]],symbol=solidcircle, symbolsize=40, scaling=constrained, axes=none, size=[100,100]), p=L):
display(Matrix(8,8,[P]));
pointplot([[0,0],seq([cos(2*Pi*k/6),sin(2*Pi*k/6)], k=0..5),[3/2,sqrt(3)/2],[2,0],[3/2,-sqrt(3)/2]], color=[red,L[28][],"LightGrey","LightGrey","LightGrey"], scaling=constrained,symbol=solidcircle, symbolsize=20, size=[600,600]);

Edit. You can get by with just one additional point above the top side of the hexagon, as  @Alfred_F  already wrote about. Regardless of the choice of its color, we get an equilateral triangle with vertices of the same color:

Look at the following 2 examples:

cat(alpha,beta); # not OK
`αβ`; # OK
cat(d,alpha); # not OK
`dα`; # OK

For your example

`#3  Q24` (`α`);

restart;
Proc := proc(m)
local xA, yA, xB, yB, xC, yC, xJ, yJ, tx, Oo, c1, r, eqBJ, eq1, sol, O;
uses plots, geometry;
_EnvHorizontalName := 'x'; _EnvVerticalName := 'y';
xJ := 5; yJ := 1; point(A, 2, 4); point(J, xJ, yJ); point(Oo, 0, 0);
r := 3; circle(c1,[Oo, r]);
eqBJ := y = m*(x - xJ) + yJ; line(BJ, eqBJ, [x, y]);
eq1 := x^2 + y^2 = r^2; sol := solve({eqBJ, eq1}, {x, y}, explicit);
xB := subs(sol[1], x); yB := subs(sol[1], y);
point(B, xB, yB); xC := subs(sol[2], x); yC := subs(sol[2], y);
point(C, xC, yC); circle(c2, [A, B, C],centername=O); line(AB, [A, B]); line(AC, [A, C]); line(BC, [B, C]);
eqBJ := y = m*(x - xJ) + yJ; line(BJ, eqBJ, [x, y]);
eq1 := x^2 + y^2 = r^2; sol := solve({eqBJ, eq1}, {x, y},explicit);
xB := subs(sol[1], x); yB := subs(sol[1], y); point(B, xB, yB);
xC := subs(sol[2], x); yC := subs(sol[2], y); point(C, xC, yC);
line(AB, [A, B]); line(AC, [A, C]);
draw([BC(color = black), c2(color = magenta), A(color = blue, symbol = solidcircle, symbolsize = 16),
B(color = red, symbol = solidcircle, symbolsize = 16), C(color = red, symbol = solidcircle, symbolsize = 16),
J(color = red, symbol = solidcircle, symbolsize = 16),c1(color=blue)],axes = normal, printtext = true);
end proc:

plots:-animate(Proc,[m],m=-0.3..1, frames=50);  # Animation

plots:-display(Proc(-0.3),Proc(0),Proc(1)); # 3 circles

It is clearly seen that these circles have only one common point  A .

Proc.mw

The formula  cos(theta) = sqrt(1 - sin(theta)^2)   is incorrect because the cosine can be negative, but the root on the right is always non-negative. The right  substitution  is  cos(theta) = sin(theta + Pi/2)  or  cos(theta) = sin(Pi/2  -  theta ) .

restart;
eqs:=[a=b,c=d,e=f]:
map(`+`, eqs, 2);
map(`-`, eqs, 2);
map(`/`, eqs, 2);
map(t->t*~2, eqs);
map(t->t^2, eqs);

I think it is a purely visual effect, because your function near  x=0  takes very small values. In fact, this function does not have a singularity at  x=0 , it is simply not defined at  0, since 0 appears in the denominator. But if you define it at x=0  by  f(0)=0 , then it will not only be continuous but also differentiable at zero. Its derivative at zero is  0 , that is, the graph touches the Ox axis. If you want to visually verify that the function is not equal to 0 near  x=0 , you can reduce the range along the vertical axis:

restart;
f:=x->piecewise(x=0, 0, sin(x)*sqrt(1 - (sin(x)/x)^2)):
f(0.1);
simplify(D(f)(x)); # Derivative of the function f 
plot(sin(x)*sqrt(1 - (sin(x)/x)^2), x = 0 .. Pi/4, y=0..0.01);

The  parametric=full  option does not work here because Maple 2018.2 does not support solving equations (inequalities) of this type. But a direct solution using  assumptions on the sign of  n  leads to success:

restart;
f := (x, n) -> 2^n/(1 + n*2^n*x^2) - 2*x^2*(2^n)^2*n/(1 + n*2^n*x^2)^2;
solve(f(x,n)>0, x, parametric=full) assuming n::integer;
Eq:=simplify(f(x,n));
solve(Eq>0, x) assuming  n>0;
solve(Eq>0, x) assuming  n=0;
solve(Eq>0, x) assuming  n<0;

Of course, this inequality can be easily solved manually, since it comes down to solving a quadratic inequality.

The minimum natural number that has the specified property is 911.
Here are some thoughts to explain the code below.
Let n=n1+n2+n3+n4+n5+n6  (the summands are in ascending order). So we have
n-n1=m1^2
n-n2=m2^2
n-n3=m3^2
n-n4=m4^2
n-n5=m5^2
n-n6=m6^2
Adding these 6 equalities we get  5*n = m1^2 + m2^2+ ... + m6^2 or
n=(m1^2 + m2^2+ ... + m6^2) /5

So we have to find 6 different squares with a minimum sum that is divisible by 5  and
(m1^2 + m2^2+ ... + m6^2) /5 > m1^2
Having tried all possible 6 of the first 30 squares, we come to a single answer:

restart;
L:=combinat:-choose([seq(i^2, i=1..30)],6):
L1:=select(t->irem(add(t),5)=0, L):
L2:=select(t->add(t)/5>t[6], L1);
L3:=[seq(add(L2[])/5-L2[][i], i=6..1,-1)];
n=add(L3);

It seems you made a typo in the expression  P(u,v) . Direct calculation shows that  P(u0,v0)<>0 . But if we replace the free term in the expression  -1  by  3  then everything is OK:

restart;
P:=(u,v)->-u^4+88*u^3*v-146*u^2*v^2+88*u*v^3-v^4+2*u^2+40*u*v+2*v^2+3;
u0 := sqrt(2)/6: v0 := -sqrt(2)/6:
P(u0,v0);
implicitdiff(P(u,v), v, u);
eval(%, [u=u0,v=v0]);
implicitdiff(P(u,v), v, u, u);
eval(%, [u=u0,v=v0]);

restart;
u:=unapply(rsolve({u(n+3)=13/12*u(n+2)-3/8*u(n+1)+1/24*u(n), u(1)=0,u(2)=1,u(3)=1}, u(n)), n);
limit(u(n), n=infinity);
sum(u(n), n=1..infinity);

There are no 4 cases here. For any value of the parameter  k , your function  x->Pre  has a single critical point  x = k/(k+1) . I do not understand why you apply the  seq command to the expression  Pre_double_prime . In this case, you simply get a sequence of operands of this expression, which is not relevant to your task.

Download Maximum_expression_new.mw

Download computation_with_dsolve_exports_new.mw

Maple easily simplifies this expression, but does not show the steps.

restart;
simplify((2+10/(3*sqrt(3)))^(1/3)+(2-10/(3*sqrt(3)))^(1/3));
#   2

But we can solve the example by hand. Let's denote this expression by  a  and cube it using the well-known formula for the cube of the sum of two numbers. Simplifying, we easily arrive at a cubic equation for a , which has a single real root a=2 :

restart;
a=(2+10/(3*sqrt(3)))^(1/3)+(2-10/(3*sqrt(3)))^(1/3):
a^3=4+3*(2+10*sqrt(3)*(1/9))^(1/3)*(2-10*sqrt(3)*(1/9))^(1/3)*a;
a^3=4+2*a;
solve(%);

If you solve this problem in Maple, then here is its full automation:

restart;
with(NumberTheory):
L:=[9,7,10,5,11]:
t1:=10*60+10: t2:=10*60+55: t3:=11*60+58:
n1:=t2-t1; n2:=t3-t2;
S1:=Divisors(n1); S2:=Divisors(n2);
S:=`intersect`(S1,S2);
select(`in`, L, S)[];

Note that  gamma  is a protected constant in Maple, not a symbol. Also, the functions  x(tau), t(tau), u(x,t) are not defined correctly.

Download non_sense_new.mw