Looks like you left out

> with(plots): 

> select(t -> (L[t] < 0), [$1..nops(L)]);

> select(t -> (L[t] < 0), [$1..nops(L)]);

The reason that

> solve(-3.75+y^1.333333333,y);

"runs amok" is that Maple is using complex numbers, and this equation is interpreted as
-3.75 + y^(1333333333/1000000000) = 0.  I suspect Maple may end up checking out 1333333333 possible solutions, corresponding to the 1333333333'th roots of 3.75, before rejecting all but one.  A more civilized answer, with only one solution, is obtained from

> solve(-3.75 + y^(4/3), y);

You could also try fsolve, which will only give you the one real solution.

I like this.  So much so that it motivated me to post an entry to my own blog, with another mathematical problem of escape from a predator.  www.mapleprimes.com/blog/robertisrael/calculusrescue

At least for the case of 3(x), I think the controlling factor is not Worksheet vs Document but the value of Typesetting:-Settings(numberfunctions).  This allows an added level of confusion, since now there is no visible clue as to what the result will be.  And a worksheet or document saved with numberfunctions=false can be opened in a session with numberfunctions=true (or vice versa), changing the results when it is re-executed.

You know that the answer should be of the form A + B*M^(2/3) for some constants A and B, but Maple doesn't.  The numerical integrator sees an expression involving a symbolic parameter and, quite correctly, refuses to handle it.  What you have to do is break the integrand up yourself.

> J:= Int(..., f = 10 .. 2670.);
  F:=  op(1, J);
  F0:= coeff(F, M^(2/3), 0);
  F1:= coeff(F, M^(2/3), 1);
  evalf(Int(F0,f = 10 .. 2670) + M^(2/3)*Int(F1, f = 10 .. 2670));

You know that the answer should be of the form A + B*M^(2/3) for some constants A and B, but Maple doesn't.  The numerical integrator sees an expression involving a symbolic parameter and, quite correctly, refuses to handle it.  What you have to do is break the integrand up yourself.

> J:= Int(..., f = 10 .. 2670.);
  F:=  op(1, J);
  F0:= coeff(F, M^(2/3), 0);
  F1:= coeff(F, M^(2/3), 1);
  evalf(Int(F0,f = 10 .. 2670) + M^(2/3)*Int(F1, f = 10 .. 2670));

For a single die, you have a "discrete uniform" distribution.  Your X is the sum of two independent discrete uniform random variables.

> with(Statistics):
   X1:= RandomVariable(DiscreteUniform(1,6));
   X2:= RandomVariable(DiscreteUniform(1,6));
   X:= X1 + X2;

Unfortunately Maple has trouble with sums of discrete random variables.  The following should work, but doesn't.

> Probability(X <= 4);

int(int(sum(1/6*(undefined*Dirac(_t7-1)+Heaviside(_t7-1)-Heaviside(_t7-6)+undefined*Dirac(_t7-6))*Dirac(_t7-k),k = -infinity .. infinity)*sum(-1/6*(-Heaviside(-1+k)+Heaviside(-6+k))*Dirac(_t8-_t7-k),k = -infinity .. infinity),_t7 = -infinity .. infinity),_t8 = -infinity .. 4)

However, you can do this:

> add(Probability(X = j), j = 1 .. 4);

            1/6

and similarly for the other questions.

For a single die, you have a "discrete uniform" distribution.  Your X is the sum of two independent discrete uniform random variables.

> with(Statistics):
   X1:= RandomVariable(DiscreteUniform(1,6));
   X2:= RandomVariable(DiscreteUniform(1,6));
   X:= X1 + X2;

Unfortunately Maple has trouble with sums of discrete random variables.  The following should work, but doesn't.

> Probability(X <= 4);

int(int(sum(1/6*(undefined*Dirac(_t7-1)+Heaviside(_t7-1)-Heaviside(_t7-6)+undefined*Dirac(_t7-6))*Dirac(_t7-k),k = -infinity .. infinity)*sum(-1/6*(-Heaviside(-1+k)+Heaviside(-6+k))*Dirac(_t8-_t7-k),k = -infinity .. infinity),_t7 = -infinity .. infinity),_t8 = -infinity .. 4)

However, you can do this:

> add(Probability(X = j), j = 1 .. 4);

            1/6

and similarly for the other questions.

You mean something like this?

> A:= Matrix(Array(0..2,0..2,`*`) mod 3):
  <<`*`|`|`|Vector[row]([$0..2])>,
   Vector[row]([`-`,`+`,`-`$3]),
   <Vector([$0..2])|Vector([`|`$3])|A>>;

                       [*    |    0    1    2]
                       [                     ]
                       [-    +    -    -    -]
                       [                     ]
                       [0    |    0    0    0]
                       [                     ]
                       [1    |    0    1    2]
                       [                     ]
                       [2    |    0    2    1]

You mean something like this?

> A:= Matrix(Array(0..2,0..2,`*`) mod 3):
  <<`*`|`|`|Vector[row]([$0..2])>,
   Vector[row]([`-`,`+`,`-`$3]),
   <Vector([$0..2])|Vector([`|`$3])|A>>;

                       [*    |    0    1    2]
                       [                     ]
                       [-    +    -    -    -]
                       [                     ]
                       [0    |    0    0    0]
                       [                     ]
                       [1    |    0    1    2]
                       [                     ]
                       [2    |    0    2    1]

Is that the end of the procedure gradPhi4?  There is no end proc.  What arguments are you calling it with?  What results do you get in Maple 11 and Maple 13?

Yes, it should give you the initial conditions, perhaps after some simplificiation.

> S:=dsolve([op(sys),op(IC)]):
   simplify(eval(S, t=0));

{Bx(0) = -exp(-1/2*x^2/sigma^2)*(-I+256*I*omega^4+32*omega^2)/(-16*omega^2+I)/(16*I*omega^2+1), By(0) = 0, Bz(0) = -3/2*I*k*x/sigma^2*exp(-1/2*x^2/sigma^2)/omega, Ux(0) = 0, Uy(0) = 0, Uz(0) = 0}
   

Almost there: needs a bit more simplification

> simplify(evalc(%));

{Bx(0) = exp(-1/2*x^2/sigma^2), By(0) = 0, Bz(0) = -3/2*I*k*x/sigma^2*exp(-1/2*x^2/sigma^2)/omega, Ux(0) = 0, Uy(0) = 0, Uz(0) = 0}

> IC;

[Bx(0) = exp(-1/2*x^2/sigma^2), By(0) = 0, Bz(0) = -3/2*I*k*x/sigma^2*exp(-1/2*x^2/sigma^2)/omega, Ux(0) = 0, Uy(0) = 0, Uz(0) = 0]

Yes, it should give you the initial conditions, perhaps after some simplificiation.

> S:=dsolve([op(sys),op(IC)]):
   simplify(eval(S, t=0));

{Bx(0) = -exp(-1/2*x^2/sigma^2)*(-I+256*I*omega^4+32*omega^2)/(-16*omega^2+I)/(16*I*omega^2+1), By(0) = 0, Bz(0) = -3/2*I*k*x/sigma^2*exp(-1/2*x^2/sigma^2)/omega, Ux(0) = 0, Uy(0) = 0, Uz(0) = 0}
   

Almost there: needs a bit more simplification

> simplify(evalc(%));

{Bx(0) = exp(-1/2*x^2/sigma^2), By(0) = 0, Bz(0) = -3/2*I*k*x/sigma^2*exp(-1/2*x^2/sigma^2)/omega, Ux(0) = 0, Uy(0) = 0, Uz(0) = 0}

> IC;

[Bx(0) = exp(-1/2*x^2/sigma^2), By(0) = 0, Bz(0) = -3/2*I*k*x/sigma^2*exp(-1/2*x^2/sigma^2)/omega, Ux(0) = 0, Uy(0) = 0, Uz(0) = 0]