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

I've tweaked your procedure a bit. Now it automatically calculates  z_min  and  z_max . I made the numerical values for the grid  as the procedure's parameter (the list  G ). Now a user can choose these values himself (without changing the procedure code) so that the graph looks beautiful :

restart:
grids:=proc(f::algebraic,x_range::`=`,y_range::`=`,G::list)
 local z_min,z_max,plot_f,xz,yz,xy;
uses plots;
z_min:=floor(minimize(f,x_range,y_range));
z_max:=ceil(maximize(f,x_range,y_range));
plot_f:=plot3d(f,x_range,y_range);
xz:=plot3d([x,op([2,1],y_range),z],x_range,z=z_min..z_max,style=line,color=blue,thickness=0,grid=[G[1],G[3]]);
yz:=plot3d([op([2,1],x_range),y,z],y_range,z=z_min..z_max,style=line,color=blue,thickness=0,grid=G[2..3]);
xy:=plot3d([x,y,z_min],x_range,y_range,style=line,color=blue,thickness=0,grid=G[1..2]);
display(plot_f,xz,yz,xy,lightmodel=none,labels=[x,y,"f(x,y)"],labeldirections=[horizontal,horizontal,vertical],axes=frame);
end proc:

# Examples
f:=x^2-y: g:=x^2-y^2:
A:=grids(f,x=-1..3,y=-2..1,[5,4,7]):
B:=grids(g,x=-1..1,y=-1..1,[5,5,5]):
plots:-display(< A | B >);

Use option AllSolutions :

restart;
Sol:=solve(sin(x) = -1/2, x, AllSolutions);
about(_B1);
eval(Sol,_B1=0), eval(Sol,_B1=1);

Do

odeplot(s, [x, abs(f(x))], x = 0 .. 5, colour = red);

Use the  assign  command:

restart;
eq_VS_m4 := 132.790562 = -0.000588*VRi + 0.995381*VRr - (1692.955600*VRi - 1747.964480*VRr)/(VRi^2 + VRr^2);
eq_VS_m5 := 0 = 0.000588*VRr + 0.995381*VRi + (1747.964480*VRi + 1692.955600*VRr)/(VRi^2 + VRr^2);
eq_VS_m6 := (VRi^2 + VRr^2) * lhs(eq_VS_m4) = simplify((VRi^2 + VRr^2) * rhs(eq_VS_m4));
eq_VS_m7 := (VRi^2 + VRr^2) * lhs(eq_VS_m5) = simplify((VRi^2 + VRr^2) * rhs(eq_VS_m5));
sol := solve({eq_VS_m6,eq_VS_m7}, {VRr,VRi});

assign(sol[3]);
sqrt(VRi^2+VRr^2);

Or by the  eval  command without assigning:

restart;
eq_VS_m4 := 132.790562 = -0.000588*VRi + 0.995381*VRr - (1692.955600*VRi - 1747.964480*VRr)/(VRi^2 + VRr^2);
eq_VS_m5 := 0 = 0.000588*VRr + 0.995381*VRi + (1747.964480*VRi + 1692.955600*VRr)/(VRi^2 + VRr^2);
eq_VS_m6 := (VRi^2 + VRr^2) * lhs(eq_VS_m4) = simplify((VRi^2 + VRr^2) * rhs(eq_VS_m4));
eq_VS_m7 := (VRi^2 + VRr^2) * lhs(eq_VS_m5) = simplify((VRi^2 + VRr^2) * rhs(eq_VS_m5));
sol := solve({eq_VS_m6,eq_VS_m7}, {VRr,VRi});

eval(sqrt(VRi^2+VRr^2), sol[3]);
eval(sqrt(VRi^2+VRr^2), sol[2]);

The method with  eval  is more flexible, because allows not to overrun the entire code to use other solutions instead  sol[3] , but to use any of the solutions found if needed.

restart;
M:=proc(M1,M2)
local L1, L2;
L1:=[seq((k-1)/(M1-1),k=1..M1)]; L2:=[seq((k-1)/(M2-1)$M1,k=1..M2)];
Matrix(M1*M2, (i,j)->SGP[j](L2[i],L1[irem(i-1,M1)+1]));
end proc:

Examples of use:

M(2,3);
M(3,2);
M(3,3);

https://www.mapleprimes.com/users/Jean-Michel%20Collard/questions

You should not use a name and the same indexed name in the same code.

Example:
 

f := a*x+b; 
f[1] := subs(x = 1, f); 
f[2] := subs(x = 2, f);

                                                    f := a x + b
                                                  f[1] := a + b
                                                     f[2] := f

Change names in your code, for example   psi[i]  with  psi1[i]  and  Omega[i]  with  Omega1[i]

We can easily find one solution to this equation in the form of a constant function  f(x) = C :

restart;
eq:=(1/24)*exp(-8)-(1/12)*exp(-5)+(1/24)*exp(-2)-(1/24)*exp(-10)+(1/24)*exp(-9)+(1/24)*exp(-1)+diff(f(x), x, x, x, x, x)+d*(diff(f(x), x))+e*f(x)+a*(diff(f(x), x, x, x, x))+b*(diff(f(x), x, x, x))+c*(diff(f(x), x, x))-1/24;
eval(eq,f(x)=C);
f(x)=solve(%, C);
odetest(%,eq);

The initial conditions for obtaining this solution are obvious. It is also obvious that  limit(C, x=infinity)=C  for C=constant