See help on the  evala  command. Also the  Algebraic  package may be useful for you.

You forgot to call the package  plots . Insert  with(plots):  at the beginning of your worksheet.

Replace in the code the dot  .  with the multiplication sign  * .

The corrected file 

test1_03_new.mw

The error you receive says that Maple does not have the  FourierSeries  package.

You can freely download the  OrthogonalExpansions  package from the Maple Application Center. This package contains a number of commands for working with Fourier series.

https://www.maplesoft.com/applications/view.aspx?SID=33406

(rhs@op)~([rootsq0]);

There will be infinitely many positive solutions too. Below we find the smallest:

restart;
eq := 61*x^2 + 1 = y^2;
[isolve(eq)];
map(t->eval([x,y],t),%);
map(t->expand(eval(t,_Z1=1)), %);
[x,y]=~select(t->t[1]>0 and t[2]>0, %)[];

                     [x = 226153980, y = 1766319049]

Edit. It is convenient to find the desired solutions using the procedure:

restart;
Eq := 61*x^2 + 1 = y^2:
Sol:=n->subsindets(eval(isolve(Eq)[4],_Z1=n),realcons,expand):

# Examples of use:
Sol(1);
Sol(2);

In fact, your example is on optimization with a parameter  a . The commands of simbolic optimization  minimize  and  maximize  can solve such examples only in very simple cases, for example  minimize(a*x+3, x=0..1);  , a little more difficult   minimize(x^2+a*x+3, x=0..1);   - it already fails. Below is a solution using a procedure  Min  that returns the result for each specific  a :

restart;
local f;
Min:=proc(a)
uses Optimization;
f:=a*x + log10(3+4.2^x +2.2^(2*x));
if a>=0 then return [evalf(eval(f,x=0)),[x=0]] else
Minimize(f, x=0..1) fi;
end proc:

# Examples:
Min(2);
Min(-2);
Min(-0.3);

You did not specify the base of the logarithm, I took it equal to  10 .

This is not an integral equation, but simply a numerical calculation of a definite integral with the parameter  t :

restart; 
rho := 1/4: mu := 1/4: T := (t-x)^rho: 
E := add(T^k/GAMMA(k*rho+mu), k = 0 .. 5);
h :=unapply(int(E*(8*x^(3/2)/(3*sqrt(Pi))+x^2+(1/2)*exp(-x)-8/3)/(t-x)^.75, x = 0 .. 1, numeric), t):

# Examples:
h(1.5);
plot(h, 1..2, labels=["t","h(t)"]);

Another way is shown below with a simple example:

restart;
Sol:=dsolve({diff(y(x),x)=x+y(x), y(0)=0}, numeric);

f:=s->eval(y(x),Sol(s));
f(1);
int('f(s)',s=0..1, numeric);

f := g->int(g(r+eps), r);

f(h);

There are no such numbers, because your equality  sqrt(a+b*sqrt(c+d*sqrt(e +f*sqrt(g)))) = h  implies that  g  is an exact square:

g=solve(sqrt(a+b*sqrt(c+d*sqrt(e +f*sqrt(g)))) = h, g);

The graph shows that the root is close to 0 :

restart;
Digits:=20:
y:=t->808.2213240*(1 - 0.63*(1993551437/1601983488 - sqrt(3)/2)^0.3)*(1 - 335345*(45188/147189 - 53/(4820*ln(2)))*335345^(131537/203808)*131537^(72271/203808)*(1 - 1/(1 + (203808*exp(-677.0138344*t))/131537)^(131537/203808))/34603964738):
plot(y, -0.1..0.1);

fsolve(y(t)=196.9594856, t=-1..1);

restart;
A:=fsolve(2*x^5-x^4-0.5*x^3+3*x^2-0.5);           
B:=fsolve(x^2-1);
sort([A,B])[];

                                  -1.,  -1.000000000,  -0.4158624444,  0.4256094086,  1.

To create an animation, it is convenient to first create a procedure that creates one frame (in the future, the procedure parameter will be an animation parameter). Below are 2 procedures. The  Cone  procedure creates the cone with the base radius  r , height  h  and base center  c . The second procedure  Plane  creates the plane  z=h-y  intersecting the axis  Oz  at a point  [0,0,h] :

restart;
Cone:=proc(r,h,c:=[0,0,0])
uses plottools, plots;
display(cone(c,r,-h, color="LightBlue"), scaling=constrained, axes=none);
end proc:
Plane:=(h,X,Y)->plot3d(h-y,x=-X..X,y=-Y..Y, color=yellow):

# Examples of use

plots:-display(<Cone(2,2)|Cone(2,3)|Cone(2,5)>);
plots:-animate(Cone,[2,h,[0,0,h]], h=1..5, frames=60);
plots:-animate(Plane,[h,1.8,1.8], h=5.7..2.5, background=Cone(2,6,[0,0,6]), frames=60);

First, in your expressions, I made some simplifying notation, for example  Cp  for  C(p[1](t))  and so on (see full list of this abbreviations  Change below). This is not necessary, but just makes it easier to visually check the simplifications made. Taking these designations into account,  A  and  B  became  A1  and  B1 . Expression  A1 has been simplified well, but  B1  is only slightly simplified, because there is no subexpression b  in it. If you want to revert to the original notation, then just do the reverse substitutions.

restart;
A := -sqrt(m(p[1](t))/m(q[1](t)))*p[2](t) - l[1]*q[1](t) + l[1]*p[1](t) + q[2](t):

B := -(-sqrt(m(p[1](t))/m(q[1](t))^3)*C(p[1](t))*m(q[1](t))^(3/2)*p[2](t)^2*m(p[1](t)) + C(p[1](t))*p[2](t)^3*m(q[1](t))^(3/2)*sqrt(m(p[1](t))/m(q[1](t))) - l[1]*p[2](t)^2*C(p[1](t))*(p[1](t) - q[1](t))*m(q[1](t))^(3/2) - sqrt(m(q[1](t)))*l[2]*(p[1](t) - q[1](t))*m(p[1](t))^(3/2) + m(p[1](t))*l[1]*p[2](t)*C(p[1](t))*(p[1](t) - q[1](t))*sqrt(m(q[1](t))) - C(q[1](t))*q[2](t)^2*sqrt(m(p[1](t)))*m(q[1](t))*(q[2](t) - 1))/(sqrt(m(p[1](t)))*m(q[1](t))^(3/2)):
indets(B);
Change:=[C(p[1](t))=Cp, C(q[1](t))=Cq, m(p[1](t))=mp, m(q[1](t))=mq, p[1](t)=p1, p[2](t)=p2, q[1](t)=q1, q[2](t)=q2]; # Simplifying substitutions
A1:=subs(Change,A);
B1:=subs(Change,B);

 
Eq1 := q1 - p1 = a;
Eq2 := sqrt(mq)*q2 - sqrt(mp)*p2 = b;

Rule1:=sqrt(a::anything/b::anything)=sqrt(a)/sqrt(b);
Rule2:=sqrt(a::anything/b::anything^3)=sqrt(a)/sqrt(b)/b;

algsubs(Eq2,normal(applyrule(Rule1, A1)));
A1:=algsubs(Eq1,%); # Simplified A1

applyrule(Rule2,B1);
B11:=normal(applyrule(Rule1,%));
B21:=numer(B11);
B22:=denom(B11);

B1:=simplify(algsubs(Eq1,B21))/B22; # Simplified B1

Subs.mw