Is the volume of the n-doughnut given by #Dⁿ(R,r) = #B^n-1(r)#S¹(R), with # being the appropriate measure.  R is radius of rotation.

There are a few problems, some syntactical, some semantical.   First let's deal with the syntax.

When you enter the code, Maple prints the error you see.  It also moves the cursor to near the offending position.  In this case it moves it just before the T in the assignment T := T + 2*sum(...).  The problem there is that the do keyword is missing. Fix that, then re-execute.  Now the cursor moves to the call to sum.  The call is syntactically incorrect, it is missing a closing parenthesis (there are also semantic issues we will deal with shortly).  Fix that and re-execute.  Now the cursor moves to the end proc.  The problem there is that you did not close the if statement.  Add an end if.  At this point we have

Simpson:=proc(f,N,a,b);
    if (N,odd)=true then
        print("N has to be even.");
    else
        N:=N:
        a:=a:    b:=b:
        h:=(b-a)/N:
        for i from 1 to N/2 do
            S := S + 4*sum(f(x2i-1)):
        od:
        for i from 1 to N/2-1 do
            T := T + 2*sum(f(x2i));
        od:
    end if;
    S:=(h*(T+S+a+b))/3;
end proc:

The procedure is now syntactically correct. You will get some warnings about undeclared locals.  Those should be declared. 

What remain are quite a few programming bugs.  Rather than me fixing them, I'll give you some suggestions:

(N,odd) is merely a sequence. You should use either type(N,odd) or N::odd. 
print("N has to be even") is okay, however, it is better to use an error statement:  error "N has to be even".  Actually, it is better to declare that the parameter N must be even in the parameter list and then do away with the internal type check.  That is, do

Simpson := proc(f, N::even, a, b)

As you will quickly find out, a := a; and b := b; do not work (they wil generate errors).  You cannot generally assign to parameters (there are exceptions, but this isn't one of them). 

The calls to sum are invalid.  Sum requires a second argument.  However, you shouldn't be using sum for definite summation.  Instead use add.  Actually, it isn't clear why you are using sum/add anyway.  You're alreading computing the sum in the loop.  The loops could be replaced with calls to add, but this is a refinement.

Is this what you want:

Matrix(3,3,f);
                                        [f(1, 1)    f(1, 2)    f(1, 3)]
                                        [                             ]
                                        [f(2, 1)    f(2, 2)    f(2, 3)]
                                        [                             ]
                                        [f(3, 1)    f(3, 2)    f(3, 3)]

Assign f appropriately
 

`&//` := proc()
local z;
    1/add(1/z, z in [_passed]);
end proc:

R_total := R1 &// (s*L1) &// (1/s/C1) &// R2;
                                            1
                     R_total := -------------------------
                                 1      1             1
                                ---- + ---- + s C1 + ----
                                 R1    s L1           R2


normal(R_total);
                                  R1 s L1 R2
                  ------------------------------------------
                                     2
                  s L1 R2 + R1 R2 + s  C1 R1 L1 R2 + L1 R1 s

collect(%,s,factor);
                                  R1 s L1 R2
                    ---------------------------------------
                     2
                    s  C1 R1 L1 R2 + L1 (R2 + R1) s + R1 R2

I don't know that Maple has many tools for directly manipulating sums.  You could always use a purely syntactical approach:

> S := Sum(kx*x[i]+ky*y[i], i=1..n);   
                                                   n
                                                 -----
                                                  \
                                            S :=   )   (kx x[i] + ky y[i])
                                                  /
                                                 -----
                                                 i = 1

> map(op(0,S), op(S));              
                                          /  n          \   /  n          \
                                          |-----        |   |-----        |
                                          | \           |   | \           |
                                          |  )   kx x[i]| + |  )   ky y[i]|
                                          | /           |   | /           |
                                          |-----        |   |-----        |
                                          \i = 1        /   \i = 1        /

> expand(%);                        
                                             /  n       \      /  n       \
                                             |-----     |      |-----     |
                                             | \        |      | \        |
                                          kx |  )   x[i]| + ky |  )   y[i]|
                                             | /        |      | /        |
                                             |-----     |      |-----     |
                                             \i = 1     /      \i = 1     /

eq := expand(cos(3*x)) = cos(3*x);
sol := solve(eq, {cos(x)}):
eval(sol[1], x=theta/3);
     theta                                      2 1/2 (1/3)
{cos(-----) = 1/2 (cos(theta) + (-1 + cos(theta) )   )
       3

                               1
     + 1/2 -----------------------------------------}
                                         2 1/2 (1/3)
           (cos(theta) + (-1 + cos(theta) )   )

Use the partition option.  That is, ApproximateInt(..., partition=5).

Fern := proc(n::posint
             , x::Array(datatype=float[8])
             , y::Array(datatype=float[8])
             , r::Array(datatype=float[8])
            )
local i, pick;
option autocompile;
    for i to n-1 do
        pick := r[i];
        if pick < 0.01 then
            x[i+1] := 0;
            y[i+1] := 0.16*y[i];
        elif pick < 0.08 then
            x[i+1] := 0.2*x[i] - 0.26*y[i];
            y[i+1] := 0.23*x[i] + 0.22*y[i] + 1.6;
        elif pick < 0.15 then
            x[i+1] := -0.15*x[i] + 0.28*y[i];
            y[i+1] := 0.26*x[i] + 0.24*y[i] + 0.44;
        else
            x[i+1] := 0.85*x[i] + 0.04*y[i];
            y[i+1] := -0.04*x[i] + 0.85*y[i] + 1.6;
        end if;
    end do;
    NULL;
end proc:

n := 10^5:
r := rtable(1..n, frandom(0..1,1), 'datatype'=float[8]):
x := Array(1..n, 'datatype'=float[8]):
y := Array(1..n, 'datatype'=float[8]):

(Fern(n,x,y,r));
plot(x,y, 'style=point', 'color=green');

Using assign prevents you from reusing the assigned variables in the set of equations.  A better approach is to symbolically solve the set of equations, then create a procedure from that solution, and plug each random combination into it.  For example,

restart;
with(combinat):

vars1 := {a, b, c, d, e, f, g, i, m}:
vars2 := [h, j, k, l, n, o, p, z]:
sol := solve({NULL
              , a+b+c+d-z = 0
              , e+f+g+h = z
              , i+j+k+l = z
              , m+n+o+p = z
              , a+e+i+m = z
              , b+f+j+n = z
              , c+g+k+o = z
              , a+f+k+p = z
              , d+g+j+m = z
             }
             , vars1):

f := unapply(map(rhs,sol) union {vars2[]}, vars2);

for cnt do
    A := randcomb(100,8);
    if nops(f(A[])) = 16 then break;
    end if;
end do:

cnt; f(A[]);

It isn't clear what you really want.  What is the expression you want to plot?  The triple integral you show fully evaluates to a number. 

I didn't see a mention of the RationalMap procedure in that paper.  Did you mean ?RationalMapImage or ?RationalMapPreImage, which exist in the ?RegularChains[ConstructibleSetTools] subpackage.

You should write this as a procedure with a single argument, the Matrix to check.  The procedure should return either true or false. Here's a straightforward implementation, with enough comments to explain what is happening:

IsDiagonallyDominant := proc(M :: Matrix)
local strict,sumrest,diag,i,j,m,n;
    strict := false;   # flag to keep track of whether at least one row is strict
    (m,n) := op(1,M);  # extract number of rows (m) and columns (n) of Matrix M
    for i to m do
        diag := abs(M[i,i]);  # Get the magnitude of the diagonal term
        # Sum the magnitudes of the other terms in the row.
        sumrest := add(abs(M[i,j]), j=1..i-1) + add(abs(M[i,j]), j=i+1..n);
        if diag < sumrest then
            return false;  # short cut; return false if row fails test
        end if;
        # Assign strict, if needed
        if not strict and sumrest < diag then
            strict := true
        end if;
    end do;
    # At this point, none of the rows fail.  
    # To be diagonally-dominant, at least one row 
    # must meet the strict equality,
    # so return the value of strict.
    return strict;
end proc:

Here's a somewhat interesting extension.  Rather than return just true or false, and fail if an element is non-numeric (or, say, Pi, which is a deficiency in the above routine), we can extend this to return a boolean expression:

IsDiagonallyDominant := proc(M :: Matrix)
local pred,strict,sumrest,diag,i,j,m,n;
    pred := true;
    strict := false;
    (m,n) := op(1,M);
    for i to m do
        diag := abs(M[i,i]);
        sumrest := add(abs(M[i,j]), j=1..i-1) + add(abs(M[i,j]), j=i+1..n);
        if diag :: realcons and sumrest :: realcons then
            if evalf(diag < sumrest) then
                return false;
            end if;
            if evalf(sumrest < diag) then
                strict := true;
            end if;
        else
            pred := pred and sumrest <= diag;
            if strict <> true then
                strict := strict or sumrest < diag;
            end if;
        end if;
    end do;
    return pred and strict;
end proc:

M := <<2|1>,<c|d>>:
IsDiagonallyDominant(M);
                                       | c | <= | d |

That is not a well-formed Maple expression.  The = and > operators have the same precedence and  are non-associative (see ?operator,precedence), so you need to add parentheses to group the components. Regardless, I don't see how false could be returned. 

How did you define the functions f[i]?  Are they really Maple procedures, for example,

f[1] := x -> x^2:

Note that E[ij] should be E[i,j]. 

myproc := proc(x0)
local i,x;
  x := x0;
  for i from 0 to 10 do
      x := x-(some stuff);
      print(x);
  end do;
end proc: