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).

You can do the same thing in Maple: Enter your expression with single forward quotes. There are certain small operations (such as arithmetic) called "automatic simplifications" that will always be done. But most operations will be suspended by the quotes. To evaluate the expression after it has been entered this way, use the ditto operator: %

'int(x, x= 0..1)';

%;

The vast majority of Maple commands use global variables (such as enviroment variables, especially Digits) to some extent. Thus they can't run multi-threaded. See ?index,threadsafe for the paltry list of exceptions.

If you are interested in a hardware floating-point answer (aka "double precision"), then the following works in parallel and should give you a processor utilization percentage in the high nineties:

Sin:= proc(x)
option threadsafe, autocompile;
     sin(x)
end proc:

Cos:= proc(x)
option threadsafe, autocompile;
     cos(x)
end proc:

SC:= Threads:-Seq([Sin(i), Cos(i)], i= 1..10^6):

Try this:

Test:= proc(j)
local x;
     plot(x^j, x= -2..2)
end proc:

for j to 3 do
     print(Test(j));
     print(`Press RETURN to continue.`);
     readline(-1);
end do:

The help at ?readline says that the -1 argument avoids the dialog box and reads from the prompt. I get a dialog box anyway.