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

What means your dual equality? The system of equations? If  only the first equation is acceptable and from the boundary conditions to retain only the first two, Maple finds the only one trivial solution. Assuming that Maple was right,  the initial boundary value problem has no solution.

dsolve({(diff(f(x), x))^2-f(x)*(diff(f(x), x, x)) = diff(f(x), x, x, x)-100*(diff(f(x), x))-0.001*(2*(diff(f(x), x))*(diff(f(x), x, x, x))-f(x)*(diff(f(x), x, x, x, x))-(diff(f(x), x, x))^2), f(0) = 0, ((D@@2)(f))(0) = 0}, f(x));

                                                                        f(x)=0