You say you want to expand K in a series in R, so you should not leave K as a constant. 

When you solve the third-order DE, you will want to just use three of the four boundary conditions.  The fourth is then used to determine the
appropriate coefficient in the series for K. So:

 series(eval(de,{F(eta,R)=add(F[k](eta)*R^k,k=0..nmax),K=add(K[k]*R^k,k=0..nmax)}),R);

 des:= [seq(coeff(%,R,k),k=0..nmax)];
 bcs:=k->(F[k](-1)=`if`(k=0,-1,0),F[k](1)=`if`(k=0,1,0),D(F[k])(-1)=0,D(F[k])(1)=0);

 for k from 1 to nmax+1 do
       F[k-1]:= unapply(rhs(dsolve({des[k], bcs(k-1)[1..3]})),eta);
       K[k-1]:= solve(bcs(k-1)[4],K[k-1])
 end do;
       

Suppose your function is f(x), 
What is the equation of the tangent line at a point (p, f(p)) on the graph?
What needs to be true for that tangent line to pass through the point you were given?

I think you mean

> p:= dsolve({deq,ci},x(t),numeric, range=0..5, output= Array([0.1*j $ j=0..50]));

There is no warning message, and as far as I can tell, the result is correct.

For example:

> plot(x*exp(-x),x=0..2);
    P:= %:
    A:= op([1,1],P);    # this should be a list of the points plotted
   maxvalue:= max(op(map(p -> p[2], A)));
    

maxvalue := .367879261799993972

> select(p -> (p[2]=maxvalue), A);

[[.999012820833333314, .367879261799993972]]

I think what you want is the output option.  For example,

> A := dsolve(desystem, numeric, output= Array([0.1*j $ j=0..10]));

The resulting A will be a 2 x 1 Array where A[1,1] gives the variable names (e.g. [ t, x(t) ]) and A[2,1] is a Matrix giving the corresponding values (in this case t in the first column and x in the second).

1) If z = g(x,y) is the equation of the tangent plane at point (1,1), you can plot z = f(x,y) and z = g(x,y) together by e.g.

> plot3d([f(x,y), g(x,y)], x = -Pi .. Pi, y = -Pi .. Pi);

If you want to tell which is which, you might give them different colours, e.g.

> plot3d([f(x,y), g(x,y)], x = -Pi .. Pi, y = -Pi .. Pi, colour=[red,blue]);

In general, if you already have two plots P1 and P2, you can put them together using display in the plots package:

> plots[display]([P1,P2]);

2) mtaylor(y^2 / x^3, [x=1,y=-1],7);

You told it what the _first_ column should be: that's what the second 1 in M(1.., 1) does.  If you want a fifth column, try

> M(1 .. , 5) := <a, b, c, d>;

Your SluiceDE contains two unassigned names, kSluice and kOrifice, which I presume should be numerical constants.  If you assign them suitable values, it ought to work.  However, it doesn't, which seems to be a bug: DEplot doesn't seem to like piecewise functions where a conditional expression depends on both independent and dependent variables.  Thus:

> DEplot(diff(y(t), t) = piecewise(y(t) > t, 1, 2), y(t), t=0..1,y=0..1);

Error, (in DEtools/DEplot/CheckDE) unable to isolate DEs for required derivatives
 

On the other hand, you can get a plot of solution curves without the direction field arrows:

> DEplot(eval(SluiceDE,{kSluice=500,kOrifice=500}), h(t), t = 0 .. 100000,
     [seq(h(0)=i,i=-9..9)], h = -10 .. 10,arrows=none);

Or you could calculate and plot the direction field "by hand".

> rs:= eval(rhs(SluiceDE), {kSluice=500,kOrifice=500,h(t)=h}):
  plot([seq(seq([[t, h],[t+4000, h+4000*evalf(rs)]],
          t= 0 .. 90000, 10000), h = -9 .. 9)], colour=black);

You assign values using :=, not =.  I mean something like this (sorry, I left out a parenthesis).

> b:= 0.3; d:= 0.2; n[0]:= 2;
   for i from 1 to 100 do  n[i]:= (1+b)*(1-d)*n[i-1] end do:
   plots[pointplot]([seq]([i,n[i]], i=0..100));

Yes, you can use a "procedure defined system" in dsolve(..., numeric).  See the help page ?dsolve,numeric,ivp..

The question is, what do you really want to do with that "argument"? 

For some purposes (in particular, simplification, or making sure something converges), assume or assuming will do what you want.  For example:

> int(x^r,x=0..1);

(-limit(x^(r+1),x = 0,right)+1)/(r+1)

> int(x^r,x=0..1) assuming r::RealRange(-1,1);

1/(r+1)

e is nothing special to Maple, just another variable.  For 2.71288... you want exp(1).  Or for the exponential function, exp(t) instead of e^t.

Also, 3.14159... is Pi, not pi.

So, for example,

> int(r^2 * exp(-2*r/a), r);

-1/4*a*(a^2+2*r*a+2*r^2)*exp(-2*r/a)

 > cos(2*Pi);

1

Well then, use that form.

> subs(x=3*t, b=5, expand(sin(x-b)));

> assume(z>2, y>z*z); is(z*z>4); is(y>z*z); is(y>4);

For me in Maple 12.01, this returns true, true and true.  I also tried Maple 6 and Maple 9.5, with the same results. 
What version of Maple are you using? 

If I understand you, A is a 17 x 6 Matrix and you would like a 51 x 2 Matrix B such that the first column of B consists of the 1st, 3rd and 5th columns of A and the second column of B consists of the 2nd, 4th and 6th columns of A.  I think the simplest way to do that is

> with(LinearAlgebra):
B := <<Column(A,1), Column(A,3), Column(A,5)>|<Column(A,2),Column(A,4),Column(A,6)>>;