A slightly smaller example:

with(LinearAlgebra):
M:=Matrix([[m,0],[0,m]]);
K:=Matrix([[k,k*sin(theta)],[cos(theta),k^2]]);
m:=.02:
theta:=.12345;
k:=10^6:
K;
M;
MatrixAdd(M,K); #Not yet fully evaluated
%; #But now it is
M+K; #Same
%; #Now full evaluation
(M+K)[1,2]; #OK
eval(M+K); #OK

You could use infnorm from the numapprox package:

ode:=diff(y(x), x, x)+(a-2*q*cos(2*x))*y(x) = 0;
sol:=dsolve(ode);
plot(eval(rhs(sol),{a=3,q=2,_C1=1,_C2=-2}),x=0..100);
M:=numapprox:-infnorm(eval(rhs(sol),{a=3,q=2,_C1=1,_C2=-2}),x=0..100,xm);
#If what you want is the maximal absolute value then you are done.
#It is taken on at x = xm, where xm is
xm;
#If you want the maximal value of the solution, then you need one more use of infnorm:
#Check if the solution is positive at x = xm
eval(rhs(sol),{a=3,q=2,_C1=1,_C2=-2,x=xm});
#In my case it was not. Otherwise you would have been done.
#Now compute infnorm of f(x) + M, which is nonnegative, then subtract M:
numapprox:-infnorm(eval(rhs(sol)+M,{a=3,q=2,_C1=1,_C2=-2}),x=0..100,xm2)-M;
#So the maximal value is taken on at x = xm2:
xm2;

You could use a loop:

N:=10:
M:=LinearAlgebra:-RandomMatrix(N,3,generator=0..1.0); #Example
K:=0.6: #Threshold
f:=(a,b)->a^2+b^3; #Example function f
S:=0:
for i from 1 to N do
  if M[i,1]<=K then S:=S+f(M[i,2],M[i,3]) end if;
end do:
   
S;

Please see my answer to your earlier question:

http://www.mapleprimes.com/questions/146778-What-Is-Dt-In-This-Expression

You are applying Dist as a function, but have not defined Dist as such. (D is differentiation).
Change to:
F := ((1/6)*z+(1/3)*z^2)/(1+(1/6)*z+(1/3)*z^2);

Dist := subs(z=t-1,F); #Change 1

RealDist := Distribution(CDF = unapply(Dist, t)); #Change 2

X:=RandomVariable(RealDist);

CommonDensity := PDF(X,t);

#Notice that Dist is negative in an interval below t = 1. Did you mean to use
F1:=unapply(piecewise(t>1,Dist,0),t);
instead of unapply(Dist, t) ?

By setting printlevel to 25 and observing the long printout it seems that what happens can be mimicked by the lines given below. I'm not saying that that is what actually happens, just that something with that effect is going on.
restart;
#printlevel:=25:
eq:=Vthl = Voh + (Vtlh - Voh) * exp(-t / t1); #Assuming that is the equation
(lhs-rhs)(eq);
expand(%,exp);
collect(%,exp);
##Now if you know that Voh is less than both of the other two, you could do
solve(eq,t);
expand(%) assuming Vtlh>Voh, Vtlh>Voh;
simplify(%) assuming Vtlh>Voh, Vtlh>Voh;

Using implicitplot:
plots:-implicitplot([x^3-x*y+y^3=7 , y=-11*(x-2)-1],x=-5..5,y=-5..5);
#Solving for y in the first equation:
res:=solve(x^3-x*y+y^3=7,y);
#Plotting results together with the line:
plot([res,-11*(x-2)-1],x=-5..5,-5..5);

If solving numerically at the end is OK then you could do something like the following.
I would leave out linalg as it is not used and shouldn't be anyway: it was replaced by LinearAlgebra way back. I have also left out RealDomain.

solve({Gl1, Gl2, Gl3, Gl4, Gl5, Gl6},{c[1],c[2],c[3],f[1],f[2],f[3]});
G78:=eval({Gl7,Gl8},%);
a[1] := 50; a[2] := 314; a[3] := -400; d[1] := a[1]; d[2] := -a[2]; d[3] := a[3]; l[1] := 314; l[2] := 293; l[3] := 150; l[4] := 415;
G78;
# A procedure which solves for psi and phi for given values of alpha and beta. Output is phi only.
#An optional third argument to p is the interval in which psi and phi are supposed to be found.
#As default interval I have (arbitrarily) chosen -Pi/2..Pi/2. You may want to change that.
p:=proc(a::realcons,b::realcons,interval::range(realcons):=-Pi/2..Pi/2) local res;
  res:=fsolve(eval(G78,{alpha=a,beta=b}),{psi=interval,phi=interval});
  if not type(res,set) then
    undefined
  else
    subs(res,phi)
  end if
end proc;
p(.1234,2.567);
p(.1234,2.567,0..2*Pi);
p(-Pi/2,Pi/2);
plot3d(p,-Pi/2..Pi/2,-Pi/2..Pi/2);

Here are some:

eq:=sqrt(a* x + b) +  sqrt(c*x + d) = m;
solve(eq,x);
indets(%[1],radical);
u:=op([1,1],%);
#Requiring u to be a square is necessary and seems to be enough:
isolve(u=q^2);
eval(eq,%);
solve(%,x);

You could use a double loop (I use numbering starting from 1):

restart;
n:=5:
M:=Matrix(n);
M[1,1]:=1:
M;
for r from 2 to n do
  for c from 2 to n do
    M[r,c]:=c*M[r-1,c]+a*M[r-1,c-1]
  end do
end do:
M;

You can use the Lyapunov function V = x^2 + y^2:

F:= [- x*(x^4 + y^4) - y,  x  - y*(x^4+y^4)];
V:=x^2+y^2;
# Vdot is
simplify(F[1]*diff(V,x)+F[2]*diff(V,y));

The rest is up to you.

There is an error: Since you have multiplied by h for each of the c's and k's you should not multiply by h at the end when you assign to y[i+1] and z[i+1].

I suggest defining y and z as vectors (and using N instead of t, which confused me):

restart;
Digits:=15:
R:= -6.478831125:
h:= 0.01:
N:= 500:
z:=Vector(N+1):
y:=Vector(N+1):
z[1]:= 0:
y[1]:= Pi/2:
for i from 1 to N do
c0:= evalf(h*R*sin(y[i])): #increment for z
k0:= evalf(h*(z[i])): #increment for y
c1:= evalf(h*R*sin(y[i]+(k0/2))):
k1:= evalf(h*(z[i]+(c0/2))):
c2:= evalf(h*R*sin(y[i]+(k1/2))):
k2:= evalf(h*(z[i]+(c1/2))):
c3:= evalf(h*R*sin(y[i]+k2)):
k3:= evalf(h*(z[i]+c2)):
z[i+1]:= evalf(z[i]+(1/6)*(c0+(2*c1)+(2*c2)+c3)):
y[i+1]:= evalf(y[i]+(1/6)*(k0+(2*k1)+(2*k2)+k3)):
end do:
#Plotting in the phase plane:
plot(y,z);
#For direct comparison (using x instead of y since y is no longer unassigned:
ode:=diff(x(t),t,t)=-6.478831125*sin(x(t));
res:=dsolve({ode,x(0)=Pi/2,D(x)(0)=0},numeric,method=classical[rk4],stepsize=.01,output=listprocedure);
X,X1:=op(subs(res,[x(t),diff(x(t),t)]));
plots:-odeplot(res,[x(t),diff(x(t),t)],0..5);
X(5);
X1(5);
y[501];
z[501];

#If you insist on a procedure, here is a minimal version:

restart;
p:=proc(inits::[realcons,realcons],h::numeric,N::posint) local i,c0,c1,c2,c3,k0,k1,k2,k3,y,z,R;
   R:= -6.478831125;
   z:=Vector(N+1):
   y:=Vector(N+1):
   y[1],z[1]:= op(inits);
   for i from 1 to N do
     c0:= evalf(h*R*sin(y[i])):
     k0:= evalf(h*(z[i])):
     c1:= evalf(h*R*sin(y[i]+(k0/2))):
     k1:= evalf(h*(z[i]+(c0/2))):
     c2:= evalf(h*R*sin(y[i]+(k1/2))):
     k2:= evalf(h*(z[i]+(c1/2))):
     c3:= evalf(h*R*sin(y[i]+k2)):
     k3:= evalf(h*(z[i]+c2)):
     z[i+1]:= evalf(z[i]+(1/6)*(c0+(2*c1)+(2*c2)+c3)):
     y[i+1]:= evalf(y[i]+(1/6)*(k0+(2*k1)+(2*k2)+k3)):
   end do;
   y,z
end proc:
#Test:  
p([Pi/2,0],0.01,500);
plot(%);

My immediate guess is that by xy you meant x*y, or if the independent variable (time) is t, then it ought to be x(t)*y(t).

If a, b, c, d are given given concrete numerical values and if an initial condition is given then you can solve numerically.

Example:

sys1 := diff(y(x), x) = a/(1-y(x))^b+c*exp(-d*x)/(1-y(x))^2; #Square brackets  omitted
sys2:=eval(sys1,{b=1,a=1,c=1,d=1});
res:=dsolve({sys2,y(0)=0},numeric);
plots:-odeplot(res,[x,y(x)],-1..0.205);

Using the formula for time corresponding to the parametric form you have here:

Int(sqrt(diff(x(t),t)^2+diff(y(t),t)^2)/sqrt(-2*g*y(t)),t=0..tau);
evalf(%);
Answer 0.5828954631