You simply need to enclose the variable name in curly braces in the call to solve: Change 

solve(ics[i], v[i])

to

solve(ics[i], {v[i]})

Here's an efficiency/stylistic point, not an error: Unless you're planning on using the values solved[i] later, there's no point in storing them in variables. It just creates uncollectable garbage. Instead, combine the solve and assign:

for i to 3 do assign(solve(ics[i], {v[i]})) end do

It is an unfortunate nuance of the syntax that any reasonable attempt to create an empty Matrix using square brackets creates a Matrix with one empty row. That is, it has one row and zero columns. Note that 

ColumnDimension([[1]]) <> ColumnDimension([[]]),

the former being 1 and the latter 0.

As far as I can tell, the only two reasonable ways to create a truly empty Matrix are

Matrix(), Matrix(0), etc.

or the more primitive

rtable(1..0, 1..0, subtype= Matrix).

These will create Matrices with zero rows and zero columns.

You wrote that it is possible to write 3*h(2)+D(h)(2) as (Psi(2)+gamma)*h(2)+D(h)(2). That is not true: Psi(2) = 1 - gamma, so Psi(2)+gamma = 1, not 3. So, the expansion of Psi(2) as 1 - gamma is irrelevant in this case. But if you did want to turn off the automatic expansion of Psi(positive rational), use 

map2(Cache:-RemovePermanent, Psi, [[1],[2]]);

`Psi/constant`:= `Psi/constant`[1], 0, `Psi/constant`[3]:

If you want to also turn off the expansion of Psi(positive rational) that is done by the expand command, use

expand(expandoff()):  expandoff(Psi):

Now execute the residue command:

residue((Psi(-z)+gamma)^2*h(z), z= 2);

     2*h(2)*Psi(3)+(D(h))(2)+2*h(2)*gamma

I see that you changed your Question, changing gamma to Eulergamma. This is wrong. Eulergamma is meaningless to Maple. The Euler constant is known to Maple simply as gamma.

Edit: Answer corrected due to a problem discovered by Markiyan below.

You were on the right track with your idea of guessing b. We just have to make that automatic and systematic. Imagine a real-valued function RSS of a real variable b. RSS returns the residual sum of squares from doing a linear regression on the logarithms to find the values of A and c. (That linear regression is done by the command Statistics:-PowerFit; you don't need to take the logarithms yourself.) Then we use Maple's numerical minimizer to find the b that minimizes RSS.

I'll assume that you can come up with lower and upper bounds for b. If not, I'll need to modify my procedure a little. In the procedure below, the first argument is your two-column Matrix, and the second argument is the range of possible b values in the form min_val..max_val. I haven't tested this procedure because I don't have any suitable data. If you send your Matrix and your b range, I'd be happy to test it.

Since t-b must be strictly positive, the upper bound for b must be strictly less than the minimum value of t. If you violate this, you'll get an error message.

RecenteredPowerFit:= proc(TW::Matrix, b_range::range(realcons))
local
     T:= TW[..,1], W:= TW[..,2],
     b:= Optimization:-NLPSolve(
          b-> Statistics:-PowerFit(T-~b, W, output= residualsumofsquares),
          b_range, method= branchandbound, objectivetarget= 0
     )[2][1],
     Res:= Statistics:-PowerFit(T-~b, W)
;
     ['A' = Res[1], ':-b' = b, 'c'= Res[2]]
end proc:

If you want the symbolic derivative to be defined at the knots, then the slopes on either side of the knot need to match exactly. Matching to 14 decimal places is not exactly. To get an exact match, use exact numbers in the input: change 2.2 to 11/5 and change 1.8 to 9/5. If you make this simple change, there'll be no undefined points for the derivative. These subleties would not make any difference if you were using numeric differentiation.

By the way, the syntax x(t):= ... is a substandard way to define a function in Maple. It can lead to trouble. The correct way is x:= t-> ....

Picking up from the end of your worksheet:

ex:= expand(%[]);
coeff(ex, w, -1)/w + coeff(ex, w, -2)/w^2;

It's a simple procedure, similiar to the Univariate that I wrote for you last night. We extract the names with indets and then subtract the parameters, which are passed as an argument.

Variables:= (e, p::{set,list,rtable}(name))->
     indets(e, And(name, Not(constant))) minus convert(p,set)
:

Usage:

Variables([(a-1)*x*y^2-b+x, x-y+x^2-c], [a,b,c]);

     {x,y}

The first argument can be anything, not necessarily a list of polynomials. The second argument could just as well be {a,b,c} or <a,b,c>. Indexed variables and parameters (such as a[0] or x[i]) are allowed without the need to declare the subscripts in the set of parameters.

We need to check that the expression has one variable and that it's a polynomial with respect to that variable.

TypeTools:-AddType(
     Univariate,
     proc(P)
     local x:= indets(P, And(name, Not(constant)));
          nops(x) = 1 and P::polynom(anything,x[])
     end proc
);

With this, we can run checks like the following, each of which returns true or false.

type(x^2+1, Univariate);
type(x*y, Univariate);

When type checks occur inside other Boolean expressions, they can be done with the `::` operator:

if (x^2+1)::Univariate then ... end if;

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)