M := proc (k::nonnegint, n::nonnegint) 

option remember;

if k = 0 or n = 0 then M(k, n) := 0 elif

k = 1 then M(k, n) := n elif

n = 1 then M(k, n) := 1 else

M(k, n) := M(k, n-1)+M(k-1, n-1)+M(k-2, n-2) end if

end proc:

Example (all values for 0<=k, n<=9):

Matrix(10, (i, j) -> M(i-1, j-1) );

 Code was edited.

To really examine your 9 graphs I removed the multiplier 100, 0 in the denominator, add the colors and change the ranges of the axes:

with(plots):

pi:=Pi:

u := sin((1/2)*pi*x)*((1-cosh(pi))*sinh((1/2)*pi*y)+sinh(pi)*cosh((1/2)*pi*y));

implicitplot([u = 10, u = 20, u = 30, u = 40, u = 50, u = 60, u = 70, u = 80, u = 90], x = -5 .. 7, y = -5 .. 7, color=[red, black, khaki, pink, yellow, blue, green, violet, gold], grid = [100, 100]);

 Edit. The set  {u = 10, u = 20, u = 30, u = 40, u = 50, u = 60, u = 70, u = 80, u = 90}  was changed to the list.

For  9 team competition in one round (round-robin tournament or all-play-all tournament)   9*8/2=36  games should be played. Here is the example of a random tournament between teams  A, B, C, E, F, G, H, K, M :

restart;

roll:=rand(0..1):

assign(seq(seq(a[i,j]=roll(),j=i+1..9),i=1..8)):

assign(seq(a[i,i]=` - `, i=1..9)):

assign(seq(seq(a[i,j]=1-a[j,i],j=1..i-1),i=2..9)):

<Vector([[``,A,B,C,E,F,G,H,K,M]])|<Matrix([[A,B,C,E,F,G,H,K,M]]),Matrix(9, (i,j)->a[i,j])>>;

Any such matrix is uniquely determined by the numbers that stand above (or below) the main diagonal. Obviously the number of the all matrices of such type is  2^36 .  If we want to write down all the possible variants in a single matrix, we get a huge matrix with  36  rows and  2^36  columns.

I think the original problem can be solved, if we go in 3-D space. We just roll up the plane in a tube (cylinder). Axis  theta goes along the circumference, and the values of the function we read from the vertical axis.

Example:

R := 1:

E := plot([cos(3*theta), 0], theta = 0 .. 2*Pi, color = [red, blue], thickness = [3, 2]):

F := plot([1, -1], theta = 0 .. 2*Pi, color = grey, filled):

plots[display](E, F, scaling = constrained);

A := plot3d([R*cos(theta), R*sin(theta), z], theta = 0 .. 2*Pi, z = -1 .. 1, style = surface, color = grey):

B := plots[spacecurve]([R*cos(theta), R*sin(theta), cos(3*theta)], theta = 0 .. 2*Pi, color = red, thickness = 4, linestyle = solid):

C := plots[spacecurve]([R*cos(theta), R*sin(theta), 0], theta = 0 .. 2*Pi, color = blue, thickness = 3, linestyle = solid):

plots[display](A, B, C, axes = normal, tickmarks = [3, 3, 3], view = [-1.7 .. 1.7, -1.7 .. 1.7, -1.7 .. 1.7], scaling = constrained);

without any equations:

plots[display](plottools[sphere]([0, 0, 0], 1));

Example:

R:=3: x0:=3: y0:=-1:  # Parameters of the circle

Circle:=plot([R*cos(t)+x0,R*sin(t)+y0, t=0..2*Pi], color=blue, thickness=2):

Data:=convert(op([1,1], plot(cos(theta), theta=0..2*Pi)), listlist):

Cosine:=plot(select(t->is((t[1]-x0)^2+(t[2]-y0)^2<=R^2), Data), color=red, thickness=3):

plots[display](Circle,Cosine, view=[min(0,x0-R)..x0+R,y0-R..y0+R], scaling=constrained);

Addition: By varying the parameters of the circle, you can get any circles. I do not understand the exact meaning of the phrase "The axis theta of the ploted function is at a radius R" 

Expression  e^x  should be coded as  exp(x)  or use  Expression  palette (in 2D Math of Standard GUI) to make it look in the usual form.  When you just type  the letter e  for Maple it just is a symbol rather than the number  e=2.7...

plots[textplot]([0, 0, "Hello"], color = green, font = [TIMES, ROMAN, 100], axes = none);

Another way - use for dummy variable any special symbol that distinguishes its special status, for example

 z := Int(f(t - _s),  _s=0..1);

You have to put the delimiters after each line (colon or semicolon), for example:

ma[1] := [-885.880598855072776, [bh = 0., g0h = 0., g1h = 0.825946224722250e-4]]:

ma[2] := [-877.957885609114328, [bh = 347.116836805625, g0h = 0., g1h = 0.164861392564e-3]]:

ma[3] := [-863.445203144006655, [bh = 0., g0h = 0., g1h = 0.787090377403668e-4]]:

aver := [bh = 121.477814540743, g0h = 0.249071340242633e-5, g1h = 0.110740782334817e-3]:

add((aver[1]-ma[j][2][1])^2, j = 1 .. 3);

You need to download commands in textual form, but not in the form of pictures, as requires a lot of time for retyping of data. Therefore, I present the solution for other data:

restart;

X:=Vector([$ 1..5]):

Y:=Vector([1.2, 3.8, 9.1, 15.7, 24.9]):

R:=unapply(Statistics[PowerFit](X, Y, t), t);

E:=add((Y[i]-R(X[i]))^2, i=1..5);  # Total Error

A:=plot(op(map(convert,[X,Y], list)), color=red, style=point, symbolsize=20):

B:=plot(R(x), x=0..5.7, color=blue, thickness=2):

plots[display](A,B);  # Visualization

Example:

from 1 to 10 do

combinat[randperm]([$ -3..3])

od;

I do not understand why you take only the first three and the last three terms of the series.

If you want to calculate  cos(70)  with 6 digits by series, then using the fact that the series for the cosine is  alternating and of Leibniz's type first find for which  n  n-th term of the series will be less than  10^(-6) .

fsolve(70^(2*n-2)/(2*n-2)! = 10^(-6), n=1..infinity);

                                 121.7465466

Therefore, it suffices to take 122 members of the series.

Two approaches. 

The first way - calculate the sum of these terms symbolically, then approximately with 6 digits:

evalf[6](add((-1)^(n-1)*70^(2*n-2)/(2*n-2)!, n=1..122));

                                 0.633319

The second method - compute all members approximately. There need high accuracy, because the first terms of the series in absolute value are very large (40 digits is enough):

Digits:=40:

add((-1)^(n-1)*evalf(70^(2*n-2)/(2*n-2)!), n=1..122):

evalf[6](%);

                                0.633319

Gauss - Seidel method is very close to the simple iterative method and differs only in that  each iteration uses just calculated components of this iteration. This method can be used for both linear and nonlinear systems, and its convergence depends on the choice of the starting point and the Jacobian of mapping, defined by the system. For details, see the paper

The basic formula for Gauss - Seidel method:

x1[n] = f1(x1[n-1] , x2[n-1] , x3[n-1])

x2[n] = f2(x1[n] , x2[n-1] , x3[n-1])

x3[n] = f3(x1[n] , x2[n] , x3[n-1])

Compare with the formula for simple iterative method:

x1[n] = f1(x1[n-1] , x2[n-1] , x3[n-1])

x2[n] = f2(x1[n-1] , x2[n-1] , x3[n-1])

x3[n] = f3(x1[n-1] , x2[n-1] , x3[n-1])

For the example above we find the initial approximation  [6, 0, -3]  of the graphs of the first two equations:

Sys:=[x1^2+x2-37=0, x1-x2^2-5=0, x1+x2+x3-3=0]:

plots[implicitplot](Sys[1..2], x1=-1..12, x2=-4..40, color=[red,blue], gridrefine=3);

From the plot clearly shows that the system has two solutions. At first, we search  the top point using the upper branch of the blue parabola:

Sys1:=[x1=sqrt(37-x2), x2=sqrt(x1-5), x3=3-x1-x2]:

f:=(x1,x2,x3)->(sqrt(37-x2), sqrt(x1-5), 3-x1-x2):

x1[0]:=6.: x2[0]:=0.: x3[0]:=-3.:

f1:=(x1,x2,x3)->f(x1,x2,x3)[1]: f2:=(x1,x2,x3)->f(x1,x2,x3)[2]: f3:=(x1,x2,x3)->f(x1,x2,x3)[3]:

for n from 1 do

x1[n]:=f1(x1[n-1], x2[n-1], x3[n-1]); # Formulas for iterations

x2[n]:=f2(x1[n], x2[n-1], x3[n-1]);

x3[n]:=f3(x1[n], x2[n], x3[n-1]);

if abs(x1[n]-x1[n-1])<10^(-5) and abs(x2[n]-x2[n-1])<10^(-5) and abs(x3[n]-x3[n-1])<10^(-5) then break fi;

od:

n, [x1[n], x2[n], x3[n]]; # The number of itarations and the solution

eval(Sys, [x1=%[2,1], x2=%[2,2], x3=%[2,3]]); # Verification

5 iterations was done for needed accuracy.

Compare the above solution with the solution by simple iterative method:

f:=(x1,x2,x3)->(sqrt(37-x2),sqrt(x1-5),3-x1-x2):

x1[0]:=6.: x2[0]:=0.: x3[0]:=-3.:

f1:=(x1,x2,x3)->f(x1,x2,x3)[1]: f2:=(x1,x2,x3)->f(x1,x2,x3)[2]: f3:=(x1,x2,x3)->f(x1,x2,x3)[3]:

for n from 1 do

x1[n]:=f1(x1[n-1],x2[n-1],x3[n-1]);

x2[n]:=f2(x1[n-1],x2[n-1],x3[n-1]);

x3[n]:=f3(x1[n-1],x2[n-1],x3[n-1]);

if abs(x1[n]-x1[n-1])<10^(-5) and abs(x2[n]-x2[n-1])<10^(-5) and abs(x3[n]-x3[n-1])<10^(-5) then break fi;

od:

n, [x1[n], x2[n], x3[n]]; 

In order to achieve the same accuracy  9 iterations required.

The second root of the above system can be found similarly, only it is necessary to take

f:=(x1,x2,x3)->(sqrt(37-x2), -sqrt(x1-5), 3-x1-x2):

Edited: header is fixed. At first I wrongly called the method of simple iteration as Newton's method.

Look, even more simple sum (similar to your sum) not simplified in either of finite or infinite range:

sum(exp(-X^2)/(1+exp(X)), X=0..n);

sum(exp(-X^2)/(1+exp(X)), X=0..infinity);

Therefore, probably the only way to compactly write  your sum is to use the inert Sum command.