If you need to write this as an inert sum (not for calculation, but for example for a presentation), then you can do it like this:

restart;
2^(k+1)*Sum(sqrt((2*r-1)*(2*s-1))*psi[n*(M+1)+s](t), s=`and`(1,(s+r)::odd) .. r-1) ;

I can not confirm your claim. In Maple 2018 we have:

0^0;
limit(x^x, x=0);
plot(x^x, x=0..2, 0..4);
                                        

PS. I just noticed that you have a very old version of Maple (20 years ago). Since then, this error has long been corrected.

You can use inert sqrt  %sqrt  or  ``(...)  construction:

ex:=sqrt(x)/2/h;
ex1:=%sqrt(ex^2);  # Or
ex2:=sqrt(``(ex^2));

#  Examples of use
simplify(ex1) assuming h>0;  # Or
value(ex1) assuming h>0;  # Or
simplify(expand(ex2)) assuming h>0;
 

A way:

restart;
assign(seq(seq(m[i,j]=x^i*y^j, j=1..5), i=1..5)):
w:=m[2,3]+m[3,2];
v:=m[1,5]-3*m[2,2]+2*m[3,4];

Output:                       
                                 
                                

Download cp.mw

Yes, we can plot the entire graph without using differential equations according to the following plan:

1. First, we plot the graphs of  f , g , h of the implicit functions, from which we find the maximum range for  t - variable, so that there are real values of  x , y , z.

2. Then, with step = 0.01, we solve the equations for these variables. For every  t  we get  2 values of  x , 2 values of  y  and 1 value of  z .

3. Thus we get 4 branches of the entire graph and plot them.

The final result:

All the codes in the attached file.

curves_4.mw

x , y, z  variables are specified by the equations f , g , h  implicitly. These equations have real solutions not for every  t = 0 .. 2*Pi. This is clearly seen from the graphs of implicit functions. For example, from the graph  f , we see that for  t>3.48...  and  t<2.2...  there are no real solutions:

plots:-implicitplot(f, t=0..2*Pi, x=-10..10, gridrefine=5, view=[0..2*Pi, -1.5..1.5]);
                         

Edit.

You can get an symbolic approximate value of this integral by expanding the integrand in a series in powers of  x :

restart;
S:=12*x^3*c[2]+6*x^2*c[1]+x^2*exp(x^3*c[2])*exp(x^2*c[1]):
P:=convert(series(S, x=0, 10), polynom);
int(P, x = 0..1);
                             

In  s:=LowerCase(s):  line of your procedure there are 2 errors at once:

1. You can not assign anything to formal parameters in the body of a procedure.

2. Strings can not act as names.
 

Here are 2 examples:
It does not work:
P:=proc(x::string)
uses StringTools;
x:=UpperCase(x);
%;
end:
P("y");

But it works:
P1:=proc(x::string)
local s;
uses StringTools;
s:=UpperCase(x);
s;
end:
P1("y");

Edit.
 

f:=x->x^2*sqrt(25-x^2):
a:=0.: b:=5.: n:=100: h:=(b-a)/n: S:=0:
for j from 0 to n-1 do
S:=S+(f(a+j*h)+f(a+(j+1)*h))/2*h;
end do:
S;

restart;
with(plots):
local gamma:
g:=2:
rr:=3:
V1:=arrow([0.5,0], [-4.5,0], width=[0.006, relative=false], head_width=[0.07, relative=false], head_length=[0.2, relative=false], view=[-4..4,-1..1]):
V2:=arrow([1.5,0],[2.5,0], width=[0.006, relative=false], head_width=[0.07, relative=false], head_length=[0.2, relative=false], view=[-4..4,-1..1]):
Seg:=plot(0, x=1/g+0.1..rr/g-0.1, color=black, thickness=7):
C:=plot([[0.55,0],[1.45,0]], style=point, color=black, symbol=circle, symbolsize=14):
L:=plot([[0,t,t=-0.04..0.04],[1,t,t=-0.04..0.04],[2,t,t=-0.04..0.04]], color=black, thickness=5):
T:=textplot([[0,-0.15,0],[1,-0.15,1],[0.55,-0.15,1/gamma],[1.45,-0.15,r/gamma],[2,-0.15,r]], font=[times,bold,16]):
display(V1,V2,Seg, C, L, T, size=[800,400], axes=none);

   Output:

restart;
y1:=-1;
y2:=eval(y, solve({y=x, x=2-y^2, x>0}));
int(1, [x=y..2-y^2, y=y1..y2]);
                                                            

x,y := 1,2;
x%+y=x+y;
                                               

Replace  sum  by  add . Now it works:

restart;                                                                  
with(plots):                                                               
g:=(x,y)->sin(Pi*(x+y))*(x^2+y^2);                                        
q:=(n,k,l,x,y)->binomial(n,k)*binomial(n-k,l)*(((1+x)/2)^k)*(((1+y)/2)^l)*(1-((1+x)/2)-((1+y)/2))^(n-k-l) ;

 K:=(n,x,y)->(((1+n)*(2+n))/4)*add(add(q(n,k,l,x,y)*int(int(q(n,k,l,t,u)*g(t,u),u=-1..-t),t=-1..1),l=0..n-k),k=0..n);  

p1:=plot3d(g(x,y),x=-1..1,y=-1..-x,color=blue):                           
p2:=plot3d(K(1,x,y),x=-1..1,y=-1..-x,color=brown):                        
p3:=plot3d(K(2,x,y),x=-1..1,y=-1..-x,color=yellow):                       
p4:=plot3d(K(5,x,y),x=-1..1,y=-1..-x,color=green):                        
p5:=plot3d(K(10,x,y),x=-1..1,y=-1..-x,color=red):                                                              
display([p1,p2,p3,p4,p5]);
 

restart;
int(eval(2*x^i*(x+n)^m*sqrt(-x^2+x), n=0), x = 0 .. 1);
int(2*x^i*(x+n)^m*sqrt(-x^2+x), x = 0 .. 1) assuming n>=1;