Note that by the substitution f(x) = diff(y(x),x), the differential equation becomes second order, and it may be solved explicitly in terms of the WeierstrassP function:

eq := diff(y(x),`$`(x, 3))-(diff(y(x), x))*(diff(y(x), x))+1 = 0;
eq := algsubs(diff(y(x),x) = f(x),eq);
dsolve(eq,f(x));

The boundary conditions then become {f(0) = 0,f(1) = 1,D(f)(0) = 0}. A solution exists if _C1 and _C2 can be properly chosen.

Unless I have misunderstood your prescription for obtaining D, it seems that the D you give does not correspond to the A you give. For the (1,1)-entry, say, of D should in my view be 2, not zero, because A^2 is the lowest A^k for which this entry is nonzero.

Anyway, the following code, I think, does what you ask for:

findD := (A::Matrix(square)) -> map(ldegree,add(x^k*A^k,k = 1..RowDimension(A) - 1)):
findD(A);

Update: I suddenly notice that it is only the non-diagonal elements that your are interested in. That, I guess, is the reason for your zero (1,1)-entry (and all the other diagonal entries as well). That, however, can be achieved by altering the above procedure to the following, say:

findD := proc(A::Matrix(square))
   local D := map(ldegree,add(x^k*A^k,k = 1..RowDimension(A) - 1));
   D - DiagonalMatrix(Diagonal(D));
end proc:
findD(A);

Have you tried, or considered, simply to convert your vector equations into sets of equations, using convert(...,set), and then taking the union of them? In order to do so, it would probably be best to write each differential vector equation in a form in which it has to equal the zero vector, i.e., to move any nonzero right-hand side to the left-hand side.

An example:

F := Vector(3,(i) -> f||i(t)):
G := Vector(3,(i) -> g||i(t)):
ode1 := map(diff,F,t) + G;
ode2 := map(diff,G,t) - F;
dsolve(
   convert(ode1,set) union convert(ode2,set),
   convert(F   ,set) union convert(G   ,set)
);

The following might do:

selectRealRows := (M::Matrix) -> Matrix(remove(has,convert(M,listlist),I)):

An example of use:

A := Matrix([
   [-1,3,2,1+I],
   [0,1,5,5],
   [1,1-I,4,3],
   [1,2,-3,2]
]);
selectRealRows(A);

First define an indexing function for a fully symmetric object:

`index/fullSym` := proc(inds::list,storage,entry)
   if _npassed = 2 then storage[sort(inds)[]]
   else storage[sort(inds)[]] := op(entry)
   end if
end proc:

Next, utilizing this indexing function, define the following procedure which finds a fully symmetric Array of arbitrary dimension and rank with zero traces (contraction of pair of indices):

findArray := proc(dim::posint,rank::posint)
   local A := Array(fullSym,(1..dim)$rank,proc() x[_passed] end proc);
   A := eval(A,solve(`if`(rank > 2,
      convert(add(A[i,i],i = 1..dim),set),   # If rank > 2
             {add(A[i,i],i = 1..dim)}        # If rank = 2
   )));
end proc:

Note that the symbol/variable name 'x' figures in front of [_passed]; it may, of course, be renamed. Note also that due to the full symmetry of the Array, it is sufficient to enforce zero trace of the first two indices only.

A simple example of use (using rank = 2 for visualization):

S := findArray(4,2);

The tilde tells Maple to apply the command argument() to each of the elements of your list. It is equivalent to

map(argument,[complexn1,complexn2,complexn3]);

Update: Kitonum beat me to the finish-line :-).

The example you provide may be done in Maple as follows:

with(DifferentialGeometry):
with(LieAlgebras):
DGsetup([t,x],M):
v := evalDG(  D_x);
w := evalDG(x*D_t);
LieBracket(v,w);

Using the while statement:

exit := false:
for i from 1 to 10 while not exit do
   print(i);
   exit := evalb(i > 4);
end do:

The boolean variable 'exit' can, of course, be made to depend on any condition(s).

A one-liner capable of handling V's with any number of elements:

seq(zip(`*`,Bits:-Split(i - 1,bits = nops(V)),V),i = 1..2^nops(V));

Perhaps the following will do:

V  := [a5,a6,a7]:
VS := []:
T := combinat:-cartprod([[1,0]$nops(V)]):
while not T[finished] do
   VS := [VS[],zip(`*`,T[nextvalue](),V)]
end do:
VS;

Note that the code works for any number of elements in V.

Update: Inspired by Kitonum's V1-map, it occurred to me that there is actually no need to use zip:

V  := [a5,a6,a7]:
VS := []:
T := combinat:-cartprod([map(x -> [x,0],V)][]):
while not T[finished] do VS := [VS[],T[nextvalue]()] end do:
VS;

NullSpace() determines a basis for the vector space that some matrix annihilates, i.e. the vector space of all vectors that some matrix sends to the zero vector. Eigenvectors() determines the set of eigenvectors of some matrix. Note that eigenvectors with eigenvalue zero lie in the null space, whereas eigenvectors with nonzero eigenvalues do not. Note also that whereas NullSpace() accepts non-square matrices, Eigenvectors() do not and should not.

Although being not quite what you are asking for, perhaps the following will do:

solve(convert(add(a||i*M[..,i],i = 1..ColumnDimension(M)),set));

It determines which linear combinations of the columns vanish. Note that linear dependency of the columns of any matrix of non-maximal rank need not be a matter of only two columns; it may depend on any number of columns.

To avoid loosing leading zeros, as you do, I would consistently implement binary numbers as lists. The Bits-package, as you use too, does exactly that. Annoyingly, though, the implementation there uses the reversed order of the bits to what I think is the standard one. Using ListTools:-Reverse, you seem to agree on this point.

The code below, although perhaps not pretty, seems to do the job, i.e., taking the two's complement of any binary number provided as a list.

twosCompl := proc(x::list({0,1}))
   local expr,compl;
   expr  := ListTools:-Reverse(x);
   compl := ListTools:-Reverse(Bits:-Split(2^nops(expr) - Bits:-Join(expr)));
   compl := [0$(nops(expr) - nops(compl)),compl[]][-nops(expr)..-1];
end proc:
# Your z4 as an example
twosCompl([1,1,0,0,0,0,1,0,0,1,0,0,0,0]);

Update: By way of the keyword parameter 'bits' in Bits:-Split, the above code may actually be simplified to the following:

twosCompl := proc(x::list({0,1}))
   local expr := ListTools:-Reverse(x);
   ListTools:-Reverse(Bits:-Split(2^nops(expr) - Bits:-Join(expr),bits = nops(x)));
end proc:
# Your z4 as an example
twosCompl([1,1,0,0,0,0,1,0,0,1,0,0,0,0]);

Building on the idea of Carl, first shift the rows, and next reverse the three columns (the bits):

B := Matrix([[0,0,1,1,0,0,1],[0,0,1,1,1,0,0],[0,1,0,1,0,1,0],[1,0,0,1,1,1,0]]):
Bnew := copy(B):
Bnew[[4,1,2,3],3..5] := Bnew[..,3..5]:   # Shift the rows
Bnew[..,[5,3]] := Bnew[..,[3,5]]:        # Reverse the bits, bit 4 staying put
Bnew;

Generally, only the last statement is returned, in your case the cube result. If you want both the square and the cube returned, you could instead do simply as follows:

Operations := (x) -> (evalf(x^2),evalf(x^3)):
Operations(2);

Note that in doing so, there is no need to define the local variables Square and Cube. But if you want these quantities explicitly named and accessible (from the outside), you could perhaps make them properties of an object (the module), doing as follows:

Operations := proc(x)
   module()
      export Square,Cube;
      Square := evalf(x^2);
      Cube   := evalf(x^3);
   end module
end proc:
Opers := Operations(2):
Opers:-Square;
Opers:-Cube;