Here's the explanation for why it appeared to work when you omitted the `>`. In that case, the call S(x) returns the list [2*x, -2*x].  Note that while it is sorted, it is not sorted numerically, but by address. That this only appears correct is because you limited the domain to positive values.  If you plotted from, say -5 to 5 you would see that the graph is odd rather than even, as it should be.

Alec's suggestions work because they delay evaluating C (and S) until x is assigned a numeric value.

By optimal do you mean the arrangement that fits the most rectangles into the larger rectangle?  Is there a restriction that requires the sides of the smaller rectangles to parallel the sides of the container? There is no general solution for this packing problem, though specific cases have been solved.  See Packing problem.

The CodeGeneration package provides a utility for translating Maple procedures to C, however, its ability is limited.

It would help if you explained what the problem is.  Is fa a driving term of a differential equation?  You would also do better to just ask about one of the problems, otherwise I might suspect that this was a homework assignment.

One way is

eq3 :=4*combine(eq1*x^3);
                           /infinity                 \
                           | -----                   |
                           |  \       (8 n + 3)     n|
                         4 |   )     x          (-1) |
                           |  /                      |
                           | -----                   |
                           \ n = 0                   /

Here's another approach

W := (+ mu[3]^3*mu[1]*sin(mu[1]*Pi)*cos(mu[2]*Pi)*sin(mu[3]*Pi)*cos(mu[4]*Pi)
      + mu[4]^3*mu[3]*cos(mu[1]*Pi)*cos(mu[2]*Pi)*sin(mu[3]*Pi)*sin(mu[4]*Pi)
      - mu[3]^3*mu[2]*cos(mu[1]*Pi)*sin(mu[2]*Pi)*sin(mu[3]*Pi)*cos(mu[4]*Pi)
      - mu[4]^3*mu[2]*cos(mu[1]*Pi)*sin(mu[2]*Pi)*cos(mu[3]*Pi)*sin(mu[4]*Pi)
     ):

Z1 := frontend(coeff,[W,mu[4],3]):
Z2 := frontend(coeff,[W,mu[3],3]):

Testzero(W - Z1*mu[4]^3 - Z2*mu[3]^3);
                                             true

(w1,k1) := selectremove(hastype, factor(Z1), `+`):
(w2,k2) := selectremove(hastype, factor(Z2), `+`):

Testzero(W - w1*k1*mu[4]^3 - w2*k2*mu[3]^3);
                                              true

Hmm, that would be 6.3e15 minutes.  Since typing isn't allowed in the input box, you'll have to hold down the up-arrow for quite a while.  It might be easier to just uncheck the Auto-Save feature.

That's a rather large expression.  You would be be better off using command line Maple and working with a Maple script.  Any easy way to do that is to remove all output (Edit -> Remove Output) then export as Maple Input (File -> Export -> Maple Input).  If you do that, you will find that you miss-entered all the V variables in the MSE expression.  That is, rather than, for example, V3, Maple is interpreting that as V*`3`. 

Following is the edited script.  Note that I've converted most of the line endings to colons (:) to suppress intermediate output.  I can execute it with command line Maple, however, I killed the computation of the diff's 'cause it was bringing my machine to its knees.

postedit: I just noticed that the MSE term is also a list (of one) rather than an expression.  I've corrected it below.

R1 := 14:
R2 := 14.2:
R3 := 20.2:
R4 := 17.5:
R5 := 18.5:
R6 := 20.7:
R7 := 22:
R8 := 18.5:
R9 := 18:
R10 := 13:
R11 := 15.5:
R12 := 15:
R13 := 12:
R14 := 15:
R15 := 21:
R16 := 30:
R17 := 25:
R18 := 37.5:
S2 := R2:
S3 := a*R2+(1-a)*(S2+T2):
S4 := a*R3+(1-a)*(S3+T3):
S5 := a*R4+(1-a)*(S4+T4):
S6 := a*R5+(1-a)*(S5+T5):
S7 := a*R6+(1-a)*(S6+T6):
S8 := a*R7+(1-a)*(S7+T7):
S9 := a*R8+(1-a)*(S8+T8):
S10 := a*R9+(1-a)*(S9+T9):
S11 := a*R10+(1-a)*(S10+T10):
S12 := a*R11+(1-a)*(S11+T11):
S13 := a*R12+(1-a)*(S12+T12):
S14 := a*R13+(1-a)*(S13+T13):
S15 := a*R14+(1-a)*(S14+T14):
S16 := a*R15+(1-a)*(S15+T15):
S17 := a*R16+(1-a)*(S16+T16):
S18 := a*R17+(1-a)*(S17+T17):
T2 := R2-R1:
T3 := b*(S3-S2)+(1-b)*T2:
T4 := b*(S4-S3)+(1-b)*T3:
T5 := b*(S5-S4)+(1-b)*T4:
T6 := b*(S6-S5)+(1-b)*T5:
T7 := b*(S7-S6)+(1-b)*T6:
T8 := b*(S8-S7)+(1-b)*T7:
T9 := b*(S9-S8)+(1-b)*T8:
T10 := b*(S10-S9)+(1-b)*T9:
T11 := b*(S11-S10)+(1-b)*T10:
T12 := b*(S12-S11)+(1-b)*T11:
T13 := b*(S13-S12)+(1-b)*T12:
T14 := b*(S14-S13)+(1-b)*T13:
T15 := b*(S15-S14)+(1-b)*T14:
T16 := b*(S16-S15)+(1-b)*T15:
T17 := b*(S17-S16)+(1-b)*T16:
T18 := b*(S18-S17)+(1-b)*T17:
V3 := S3+T3:
V4 := S4+T4:
V5 := S5+T5:
V6 := S6+T6:
V7 := S7+T7:
V8 := S8+T8:
V9 := S9+T9:
V10 := S10+T10:
V11 := S11+T11:
V12 := S12+T12:
V13 := S13+T13:
V14 := S14+T14:
V15 := S15+T15:
V16 := S16+T16:
V17 := S17+T17:
V18 := S18+T18:
MSE := ((R3-V3)^2
        +(R4-V4)^2
        +(R5-V5)^2
        + (R6-V6)^2
        + (R7-V7)^2
        + (R8-V8)^2
        + (R9-V9)^2
        + (R10-V10)^2
        + (R11-V11)^2
        + (R12-V12)^2
        + (R13-V13)^2
        + (R14-V14)^2
        + (R15-V15)^2
        + (R16-V16)^2
        + (R17-V17)^2
        + (R18-V18)^2)*(1/16):
diff(MSE, a);
diff(MSE, b);

A few comments.  First, you should use Array rather than array, the latter is older structure.  I'd probably use a table rather than an Array, just because ...

To generate a point plot you create a listlist structure of the form [[x1,y1],[x2,y2],...].  In this case, however, because the x are merely indices, it is easier to use Statistics:-LineChart, which only asks for a list (or Vector) or points:

L := table():
x:=0.3:
for i to 30 do
    x := 4*x*(1-x);
    L[i] := x
end do:
L := convert(L,'list'):
Statistics:-LineChart(L);

You mean, I presume, for a relative path.  Isn't that the case?  Or maybe that was what you were saying, somewhat obliquely (the "should" confused me).

I believe so. I'm surprised that fact is not mentioned in the help page, I'll submit that as a bug.

alias(seq(x=x(t), x in {a,v,p})):
deqs := eval({a = diff(v,t)
              , v = diff(p,t)
             }, a = cos(t)):
ics := eval({p = 3
             , v = 0
            }, t = 0
           ):
dsolve(deqs union ics);
                         {p = -cos(t) + 4, v = sin(t)}

Sometimes it is useful to apply a single operation to multiple operands of an expression.  The applyop procedure doesn't work here; it can apply an operation to multiple operands, but separately.  For example

applyop(f, {1,-1}, a+b+c+d);
                              f(a) + b + c + f(d)
 

I'm thinking of an applyops procedure that does

applyops(f, {1,-1}, a+b+c+d);
                               f(a+d) + b + c

Similarly, it would be useful to be able to apply an operation to just the numerator or denominator of an expression.

The integral in question appears to be

A := Int(Int(q*exp((1-b)^2*qz^2/lambda2)*exp((1-b)^2*q^2/lambda1)/Pi/(1+b^2*q^2+b^2*qz^2)
        , q=-infinity..infinity)
    , qz = 0..infinity);

The notation used, however, is not entirely clear.  That is, it uses exp^some_val; I assume that should be exp(some_val). Evaluating that directly in Maple gives 0 because the integrand is odd in q.  Should the limits of the inner integrand be symmetric?

Note that Maple uses Pi, not pi, to stand for the symbolic constant 3.14159...