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.

L:=[]: x[0]:=0:

for i from 1 to 20 do

x[i]:=fsolve(KummerM(1/2-1/4*sqrt(2*Nu),1,sqrt(2*Nu)),Nu=x[i-1]+0.1..infinity);

L:=[op(L),x[i]]:

od:

L;

I think that for all symbolic computation of multiple integrals basic  int  command is called.
At first I tried to calculate this integral using  int  command as iterated integral in 2 ways:

int(2*x+2*y+2*z, [z=-sqrt(1-(1/4)*x^2-(1/4)*y^2)+1 .. sqrt(1-(1/4)*x^2-(1/4)*y^2)+1, y=-sqrt(-x^2+4) .. sqrt(-x^2+4), x=-2..2]);    # The first way

int(int(int(2*x+2*y+2*z, z=-sqrt(1-(1/4)*x^2-(1/4)*y^2)+1 .. sqrt(1-(1/4)*x^2-(1/4)*y^2)+1), y=-sqrt(-x^2+4) .. sqrt(-x^2+4)), x=-2..2);    # The second way of the same

In the first way, the result is right, in the second -  incorrect. Probably the difference is that in the first way, the restriction   x=-2..2  is taken into account from the outset.

To find out at which step in the second method error occurs, I did all the steps of calculating:

restart;

Int1 := expand(int(2*x+2*y+2*z, z = -sqrt(1-(1/4)*x^2-(1/4)*y^2)+1 .. sqrt(1-(1/4)*x^2-(1/4)*y^2)+1));

Int2 := int(Int1, y = -sqrt(-x^2+4) .. sqrt(-x^2+4));  # Incorrect result  

Int3 := int(Int1, y = -sqrt(-x^2+4) .. sqrt(-x^2+4))  assuming x >= -2, x <= 2;  # Correct result  

int(Int2, x = -2 .. 2);

int(Int3, x = -2 .. 2);

Obviously integral  Int2  calculated incorrectly.

Conclusion: more reliably to calculate multiple integrals symbolically by syntax

int(f(x,y,z), [z=..., y=..., x=...])  if iterated integral are calculated in that order. 

restart;

X:=l1*cos(theta1)+l2*cos(theta1+theta2)=x:

Y:=l1*sin(theta1)+l2*sin(theta1+theta2)=y:

S:=solve({X, Y}, {theta1,theta2}):

simplify([allvalues(eval(S, {l1=5, l2=5, x=5, y=3}))]);  # Symbolic solution

evalf(%);  # Numeric solution (angles in radians)

Your equation is singular at the  r=0  of the factor  1/r. The equation was solved numerically for specific values of the parameters and replacement  0  by the small number  0.0001 .

restart;

n := 1: Nu := 4:

eq := diff(T(r), r, r)+(diff(T(r), r))/r+(3*n+1)*Nu*(1-r^((n+1)/n))*T(r)/(n+1) = 0;

bc := T(1) = 0, D(T)(0.0001) = 1;

Sol := dsolve({eq, bc}, T(r), numeric, maxmesh = 2056, abserr = 10^(-3));

plots[odeplot](Sol, [r, T(r)], r = 0.0001 .. 1);

You can not use the name  I[b], because  I  is the protected constant (imaginary unit). So I[b] replaced with M .

Example:

M:=<1,2; 2,3>;  Q:=<-2,3; 1,5>;

Q^(%T).M.Q;

IsMatrixShape(%,symmetric);

Maybe  OrthogonalExpansion  package help you. It  can be downloaded from 

http://www.maplesoft.com/applications/view.aspx?SID=33406

menu Tools -> Options -> Display -> Uncheck "Automatically display legends in 2-D plots"

simplify(conjugate(exp(I*x)))  assuming x::real;

Your example:

X:=[1,2]:

Y:=[[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29], [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]]:

Z:=[[20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49], [21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50]]:

plots[surfdata]([seq([seq([X[i],Y[i,j],Z[i,j]],j=1..nops(Y[i]))], i=1..nops(X))], axes=normal, view=[-0.5..2.5,-0.5..30.5,-0.5..50.5], orientation=[-40,70], labels=[x,y,z]);

decToBin:=proc(n::integer)

if n<0 then return NULL fi;

decToBin(0):=0;

if n>0 then decToBin(n):=parse(cat(decToBin(iquo(n,2)), irem(n,2)))  fi;

end proc:

Examples of use:

decToBin(10);  decToBin(163);  decToBin(-1);  decToBin(A);

restart;

f:=n->2*(1-n)/3:

k:=0: V:=Vector():

for i from -5 to 5 do

if type(f(i), integer) then k:=k+1;  V(k):=f(i)  fi;

od:

V;

I also do not understand why you introduce vectors, if at the beginning of your worksheet the partition is defined through lists. So I replaced vectors with the lists.

The following code on each call returns a random partition of  the number  n :

Partition := proc (n::posint)

local i, Integer, Sum, Partition;

Sum := 0; Partition := [ ];

for i to n-1 do

if Sum < n then Integer := RandomTools[Generate](integer(range = 1 .. n));

if Sum+Integer <= n then Partition := [op(Partition), Integer];

Sum := Sum+Integer end if end if;

end do;

`if`(Sum < n, sort([op(Partition), n-Sum]), sort(Partition));

end proc:

Example of use (the procedure was called 100 times):

seq(Partition(5), i = 1 .. 100);

Download Graph.mw