S:= A*cos(m*x)*sin(n*z) + B*cos(k*z)*sin(q*x) + sin(s*z) + C*cos(p*z)+E*sin(t*z)^2: 

L:=[]:

for i in [op(S)] do

if (nops(i)=1 and (is(op(0,i)=sin) or is(op(0,i)=cos))) or

(nops(i)=2 and (is(op([-1,0],i)=sin) or is(op([-1,0],i)=cos)) and

(type(op(1,i),constant) or type(op(1,i),symbol))) then

if nops(i)=1 then L:=[op(L), [1,i]] else L:=[op(L), [op(1, i),i]]: fi:

fi:

od: 

L;

                                          [[1, sin(s*z)], [C, C*cos(p*z)]]

A:=Array([1,1,1,0,0,3,2,4]):
n:=ArrayNumElems(A):
mean:=add(A[i], i=1..n)/n;
variance:=add((A[i]-mean)^2, i=1..n)/n;

                     mean := 3/2
                   variance := 7/4

Good question!

We will consider the integers from 0 to 10^23 -1. We assume that each of them is written by 23 digits (from 0 to 9). If the integer is less than 10 ^ 22, then write forward the required number of zeros.
Procedure   Total(N)  finds the total number of numbers in the specified range  0 .. 10^23 -1 , whose sum of the digits is equal to  N .

Total:=proc(N::nonnegint)

local L, i, j, k, l, m, n, a, b, M, s, S, T;

if N>207 then return 0 else

L:=[seq([i,0], i=0..54)]:

for i from 0 to 9 do

for j from 0 to 9 do

for k from 0 to 9 do

for l from 0 to 9 do

for m from 0 to 9 do

for n from 0 to 9 do

a:=i+j+k+l+m+n:

L:=subsop(a+1=[a,L[a+1,2]+1], L):

od: od: od: od: od: od:

M:=[seq([i,0], i=0..45)]:

for i from 0 to 9 do

for j from 0 to 9 do

for k from 0 to 9 do

for l from 0 to 9 do

for m from 0 to 9 do

b:=i+j+k+l+m:

M:=subsop(b+1=[b,M[b+1,2]+1], M):

od: od: od: od: od:

s:=0:

for i from 0 to 54 do

for j from 0 to 54 do

for k from 0 to 54 do

for m from 0 to 45 do

if i+j+k+m=N then s:=s+1; S[s]:=[i, j, k, m]; fi;

od: od: od: od:

S:=[seq(S[t], t=1..s)]:

T:=0;

for i in S do

T:=T+L[i[1]+1,2]*L[i[2]+1,2]*L[i[3]+1,2]*M[i[4]+1,2];

od;

fi;

T;

end proc:

Solution of the initial problem:

add(Total(N), N=10..99);

add(Total(N), N=100..207);

38645951372988192979800
61354048627011778971400

I guess you consider the function of two variables as a scalar field in the plane and the stream lines - the curves in the plane, which at each point has a tangent vector, which coincides with the gradient. If so, then your problem can be solved as follows (psi=3):

psi:=(x,y)->2*y-exp(-x):

C:=3: X:=[seq(i, i=0..1, 0.1)]:

Y:=[seq(solve(psi(X[i], y)=C), i=1..11)]:

V:=[seq(dsolve({diff(x(t), t)=subs(x=x(t), y=y(t), diff(psi(x,y), x)), diff(y(t), t)=subs(x=x(t), y=y(t), diff(psi(x,y), y)), x(0)=X[i], y(0)=Y[i]}), i=1..11)]:

Sol:=evalf(subs(_Z1=0, V));

plot([seq([rhs(Sol[i,1]), rhs(Sol[i,2]), t=0..3], i=1..11)], color=red, thickness=2, labels=[x, y], title="Stream lines for psi =3,  x = 0 .. 1", titlefont=[TIMES, ROMAN, 18], view=[0..2, 0..5+C] ); 

Maple does this automatically by  combine  command. But the inverse convertation of the sum to the product - it is a problem:

combine(cos(3*x)*cos(x));
combine(cos(3*x)*cos(x) + sin(3*x)*sin(x));

              1/2*cos(2*x)+1/2*cos(4*x)
                            cos(2*x)

L:=[f(-6)=2400, f(-4)=432, f(-3)=120, f(-2)=16, f(-1)=0, f(0)=0, f(1)=-8, f(2)=0]:
M:=[seq([op(lhs(L[i])), rhs(L[i])], i=1..nops(L))];
CurveFitting[PolynomialInterpolation](M, x);

           M := [[-6,2400], [-4,432], [-3,120], [-2,16], [-1,0], [0,0], [1,-8],[2,0]] 

                                               2*x^4-6*x^2-4*x

I noted  the entries in lowercase instead of uppercase. Due to the associative law it can be counted easily:

A:=Matrix([[a, b], [c, d]]):
B:=Matrix([[0], [1]]):
(A^6).B;

To solve these problems, use Mathematica rather than Maple:

Reduce[{0 < 2*f[1], 5 < f[1] + f[2], 8 < f[1] + f[3], 
 f[1] + f[4] == 7, 0 < 2*f[2], f[2] + f[3] == 6, 6 < f[2] + f[4], 
 0 < 2*f[3], 11/2 < f[3] + f[4], 0 < 2*f[4]}, {f[1], f[2], f[3], 
 f[4]}]

                                         False 

The solution will be correct, if we write

solve({ a + b = 1, a > 0, b > 1/2}, {a,b,c})  assuming  c>0;

               {a = -b+1, c = c, 1/2 < b, b < 1}

Numbers  (-1) ^ (1/2 ) ,  (-1) ^ (3/2)  and so on - are some complex numbers. If you want them to have remained in that state, and not calculated to form  a+I*b (I - complex unit), you can write like this

sum(-y2*y3+8*(-1)^``((k-1)*(1/2))*sinh(y2*Pi*k)*sin(k*Pi*y3)/(Pi^3*k^3*cosh((1/2)*k*Pi)), k = 1 .. 10);

{p2=0.00251}, {p2=1.59454}: 

assign(op(%[1])):

p2;

           0.00251

You can plot the indefinite integral:

M:=10^8:
plot(int(x^2*exp(Pi*x/10^12*M), x), x = 0 .. 1000);

To better compare the graphics,  scaling=constrained option is used and the range on vertical axis is limited.

V:=(2*x^2+2*y+gamma)*y^3:

a0:=subs(x=0.9,gamma=0,V):

a1:=subs(x=0.9,gamma=2,V):

a2:=subs(x=0.9,gamma=4,V):

b0:=subs(x=0.1,gamma=0,V):

b1:=subs(x=0.1,gamma=2,V):

b2:=subs(x=0.1,gamma=4,V):

A:=plot([a0,a1,a2], y = 0 .. 1, 0..5, linestyle = [dot, dash,dot], color = [red, blue,green]): 

B:=plot([b0,b1,b2], y = 0 .. 2, 0..5, linestyle = [dot, dash,dot], color = [red, blue,green]):

plots[display](< display(A) | display(B) >, scaling=constrained);

a1:=plot(x^2,x=1..2): a2:=plot(x^2,x=3..4): a3:=plot(x^2,x=4..5):

b1:=plot(x^3,x=1..2): b2:=plot(x^4,x=1..2): b3:=plot(x^5,x=1..2):

a:=a1,a2,a3:  b:=b1,b2,b3:

A:=Vector([a]):

B:=Vector([b]):

with(plots):

display(< A|B >);

display(LinearAlgebra[Transpose](< A|B >));

P:=proc(L::list)

local M, F, Max, S, i;

M:=convert(L, set);

F:=[seq([M[i], ListTools[Occurrences](M[i], L)], i=1..nops(M))];

print(F);

Max:=max(seq(F[i, 2], i=1..nops(M)));

S:=[];

for i in F do

if i[2]=Max then S:=[op(S), i[1]]; fi;

od;

op(S);

end proc: 

Example:

L:=[3,3,3,3,4,4,4,3,3,3,3,3,3,3,3,2,4,3,1]:

P(L);

      [[1, 1], [2, 1], [3, 13], [4, 4]]

                       3