First off, I have to say that I don't know what it means when an integral's variable of integration is also used in the limits of integration, especially when the you take the derivative of that integral with respect to that variable. I don't know what it means mathematically, nor do I know how Maple interprets it (the derivative part is especially confusing to me). However, Maple seems to not be bothered by it, so I'll work with it to give you an Answer.

We need to do the inner integral also with capital I so that they are done numerically. The differentiation will still be done symbolically before the integration. We will pass the epsilon option to all the integrals to reduce the requested accuracy. But this has to be done after the differentiation because extra integral arguments cause confusion for diff. That is what the subsindets command is doing. We cannot use method _d01ajc. I'll use the default method.

H3p:= proc(ms::algebraic, b::numeric, eps::positive:= 0.5*10^(1-Digits))
local x:= indets(ms, And(name, Not(constant))), D;
     if nops(x)>1 then
          error "1st argument must be an expression of at most 1 variable"
     end if;
     x:= x[];
     D:= diff(ms,x);
     evalf(
          subsindets(
               Int(
                    diff(
                         D*Int((-b*x+1)*Int(D^2, x= 0..x), x= 1..x),
                         x
                    )*ms,
                    x= 0..1
               ),
               specfunc(anything, Int),
               J-> Int(op(J), epsilon= eps)
          )
     )
end proc;

This procedure will give you five digits of precision in 13 minutes and ten digits in 25 minutes.

CodeTools:-Usage(H3p(ms, 0.7, .5e-5));

     memory used=58.05GiB, alloc change=488.00MiB, cpu time=12.60m, real time=12.45m, gc time=2.01m
     10.7274995619619 + 2.56038670633894*10^(-8)* I

It is simply

Int(t^2*exp(-s*t), t= 0..3);

which can be done by two easy integration by parts. (For hand computation, I use the "tabular integration" method.)

You simply need some square brackets to enclose the list of plots:

plot([plo[1], plo[2]]);

You can plot all the solutions together on the same axes with

plot(convert(plo, list));

Use a loop to find the min, like this:

Min:= infinity:
for k to 3 do
     g[k]:= DirectSearch:-DataFit(m3, model[k]);
     if g[k][1] < Min then Min:= g[k][1]; Min_k:= k end if
end do;

plots:-display(
     plot(m3, style= point, symbol= diamond, symbolsize= 9),
     plot(eval(model[Min_k], g[Min_k][2]), x= 0..27, color = black)
);

The bug is in the MTM version of dsolve. So, use the regular dsolve. You'll need to express the initial conditions as equations, not just values, so change your defining of inits to 

inits:= seq(y[k](0) = Y0[k,1], k= 1..n);

And change the dsolve command to 

Sol:= :-dsolve([odesys, inits], [seq(y[k](t), k= 1..n)], convert_to_exact= false);

Note carefully the colon-dash in front of dsolve: :-dsolve. That forces it to use the default, global dsolve rather than MTM:-dsolve.

Be prepared to either wait days or weeks (maybe more) for a solution or run out of memory first. And I'm not guaranteeing that it will ever find a solution, even if it does have enough memory. 

Since you are integrating very "choppy" functions (lots of Heaviside) over hyperrectangles in three dimensions and you only require low accuracy, try the Cuba library, new to Maple 18. See ?evalf,int,cuba. There are four Cuba methods. Each has many options, the adjustment of which will affect the speed.

First, please indent your code something like the following. It makes it much easier for me to read and debug.

f := proc(z::nonnegint, n::posint)
local x, i, k, j;
     for i to n do
          if i = n then x[i] := trunc(z/10^(n-1))
          else x[i] := trunc((`mod`(z, 10^i))/10^(i-1))
          end if
     end do;
     printf("The number %d has the digits:", z);
     x[n+1] := 0;
     for k from n by -1 to 1 do
          print(x[k]);
          x[n+1] := x[n+1]+x[k]
     end do;
     printf("The checksum is:");
     print(x[n+1]);
     for j to n do
          if x[j] <> 0 and `mod`(z, x[j]) = 0 then printf("The number is divisible by: %d\n", x[j])
          else printf("The number is not divisible by: %d\n", x[j])
          end if
     end do;
     for j to n do
          if x[j] <> 0 and `mod`(x[n+1], x[j]) = 0 then printf("The sum is divisible by: %d\n", x[j])
          else printf("The sum is not divisible by: %d\n", x[j])
          end if
     end do
end proc;

The code above has the bug fix: Before doing mod, check if the divisor is 0. That's why I have x[j] <> 0 and .... If you are using 2D input in Maple 18 or earlier, the above fix will not work. Either switch to 1D input or let me know and I will give you another fix.

Note that I declared types for your parameters: z::nonnegint and n::posint.

There are other stylistic changes that I would make: Look up the help for the integer division commands iquo and irem (?iquo)

What do want to do with the expression? It has 44,784 terms (as determined by nops). You can see any individual term, for example the 1000th term, by

op(1000, Eq1);

1. What do you mean by "eliminate"? Do you mean to set to 0 or 1? You can find all derivatives taken with respect to x in an expression ex by

indets(convert(ex, diff), 'diff'(anything, identical(x)));

You can find all derivatives of any abstract function of the three variables x, y, z, in that order, taken with respect to x in an expression ex by

indets(convert(ex, diff), 'diff'(anyfunc(identical~([x,y,z])[]), identical(x)));

2. The gradient of a vector-valued function of vector arguments is a matrix called the Jacobian. You can get it with command VectorCalculus:-Jacobian.

Here's a binary infix operator that meets your needs. I wrote it to work for Vectors or lists. In order to work properly with sort, the operator must return true when the vectors are equal.

`&<`:= proc(x::~Vector, y::~Vector)
local n:= op(1,x), i;
     for i to n while x[i] = y[i] do end do;
     i > n or i::even and x[i] > y[i] or i::odd and x[i] < y[i]
end proc;

You can use it like ordinary <:

[1,0,0,1] &< [1,1,0,1];

To use it to sort a list of vectors, pass it to sort with the quotes:

sort(L, `&<`);

I coded this under the assumption that the vectors would be the same length. If you need it to work on vectors of different lengths, let me know.

Your post is very difficult to follow. I think that your main problem is your failure to use a multiplication sign after the C. With that change, it is very easy to get nonlinear fits, as in the following worksheet. It looks like the exponential fit is better.

Download NonlinarFit.mw

You have many problems in the code. Some are errors, and some are just sloppiness. As Georgios already mentioned, Restart should be restart, and gamma cannot be global.

Some other problems:

The proc line should not end with a semicolon.

The local variables should be declared local.

Your algebraic expressions are very bloated with unnecessary parentheses.

The local value sigma is assigned but never used.

Here's my slightly improved code:

restart;
global a:=0.0065, R0:=1.225, T0:=273.16, R:=287.04;
Spd:= proc(H, T1, IAS)
local gamma:= 1.4, T, rho, sigma, EAS, TAS, SpdSnd;
     T:= T1+T0;
     rho:= R0*(1-(a/T*1000)^4.26);
     sigma:= T0/T*(1-(a/T0*H)^5.26); #Warning: not used!
     EAS:= sqrt(rho*IAS^2/R0);
     TAS:= EAS*sqrt(R0/rho);
     SpdSnd:= sqrt(gamma*R*T);
     TAS/SpdSnd;
end proc;

Another method: Use a sparse table:

 T := table(sparse, ["green" = "gruen", "red" = "rot", "blue" = "blau"]);

Then unassigned table entries will equal 0.

I've compared the two methods timewise, and I find no significant difference when "hits" are unlikely. (The time difference is within the amount I usually atrribute to random variation.)

restart:
N:= 2^22:
for k to 3 do L||k:= RandomTools:-Generate(list(integer, N)) end do:
ct:= 0:
st:= time():
T:= table(sparse, L1=~L2):
for k in L3 do if T[k] <> 0 then ct:= ct+1 end if end do:
time()-st;
ct;

     25.641
     14

restart:
N:= 2^22:
for k to 3 do L||k:= RandomTools:-Generate(list(integer, N)) end do:
ct:= 0:
st:= time():
T:= table(L1=~L2):
for k in L3 do if assigned(T[k]) then ct:= ct+1 end if end do:
time()-st;
ct;

     25.031
     14

I also did a test where hits were likely, and in that case using assigned was significantly faster. The only difference in the code is that integer becomes integer(range= 1..N). This makes the probability of a hit 1 - exp(-1) ~ 63%.

That should be dsolve({sist, cond[i][]}, ...) rather than dsolve({sist, cond[i]}, ...). Notice the extra empty square brackets. They are to convert the list cond[i] into a sequence.

There is already a command that is fairly close to what you want: numtheory:-factorset. You just need to multiply the elements:

rad:= (n::integer)-> `*`(numtheory:-factorset(n)[]);

Answering your followup question, note that the product of a list or set L is `*`(L[]), not `*`(L).