Maple doesn't solve such equations in matrix notations. If the matrices  B  and  C  are known, the equation can be easily reduced to a system of polynomial equations that Maple can solve numerically.

An example:

A:=Matrix(2, symbol=x);

B:=Matrix(2, [-1,2,3,4]):

C:=Matrix(2, [-1,0,4,-5]):

Sys:={seq(seq((A^2+B.A+C)[i,j], j=1..2), i=1..2)};

evalf(allvalues(solve(Sys)));

The equation has no solution. 

Here are my arguments:

1. Your equation does not contain the unknown function  y(x) (only its derivatives) and your initial conditions  do not contain any conditions on  y(x)  itself. Therefore, if the equation has a solution  y(x) , then any function  y(x)+constant would also be a solution. 

2. I added the condition  y(0)=0  and dropped the condition  D@@2)(y)(0) = 0 .  We got  a unique solution. But the plot of the second derivative clearly shows that  (D@@2)(y)(0) <> 0 :

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

cond := y(0)=0, (D(y))(0) = 0, (D(y))(1) = 1;

de := dsolve({cond, eq}, y(x), numeric);

plots[odeplot](de, [x, diff(y(x),x,x)], x=0..1);

Addition. Without any plot:

de(0):

                      [x = 0., y(x) = 0., diff(y(x), x) = 0., diff(y(x), x, x) = 1.39362859931605]

Maybe Lyapunov function will help you.

An example for lists:

M:=[seq(i^2, i=1..20)]:  # list 1

N:=[seq(i^3, i=1..20)]:  # list 2

plot(M,N);  # or

plot(M,N, style=point);

The same example for vectors:

M:=Vector(20, x->x^2):

N:=Vector(20, x->x^3):

plot(M,N);  # or

plot(M,N, style=point);

All the types of arrows you can see in the help on  DEtools[DEplot].

For example:

DETools[phaseportrait]([D(x)(t)=-0.4*x(t)+(0.5+4*x(t))*y(t),D(y)(t)=0.4*x(t)-(4.5+4*x(t))*y(t)],[x(t),y(t)], t=0..100, [[x(0)=1,y(0)=0.0]], x=0..1,y=0..0.1, stepsize=0.01, scaling=UNCONSTRAINED, linecolour=BLACK, dirgrid=[17,17], linestyle=1, arrows=mediumfill, axes=BOXED);

If we assume that the parameters  r  and  theta  are real numbers,  evalc  command calculates the initial expression to the complex number written in algebraic form:

sqrt(1-r^2/exp(2*I/theta));

simplify(evalc(%));

 Edited  (correction):

simplify(evalc(sqrt(1-r^2*exp(2*I*theta))));

This is not the error but a warning:

restart:

with(plots):

Digits:=35:

ini1:=D(x)(0)=0, x(0)=1:

dsys:=diff(x(t),t,t)+(x(t)-2)*diff(x(t),t)+5*x(t)=0;

dsol1 :=dsolve({dsys,ini1}, numeric, abserr=1e-9, relerr=1e-8, maxfun=0);

plots:-odeplot(dsol1, [x(t), diff(x(t),t)], t=0..2, axes=boxed, color=black, linestyle=1, tickmarks=[6, 6], axes=boxed, titlefont=[SYMBOL,12]);   # The plot  x'(t) against  x(t)

The issue is the appearance of multiple roots. See

solve((x+y+z)^2=0, z);

                                   -x-y,  -x-y

A workaround:

plot3d(solve((x+y+z)^2=0, z, dropmultiplicity = true), x=-100..100, y=-100..100);

A:=(cos(sqrt(2)*t)-cos(omega*t))/(omega^2-2);

``(1/denom(A))*numer(A);

For better visibility animation done in 2D. It is easy to rewrite in 3D by using circular symmetry about Ox-axis.

A := plot([-1, 1], x = -10 .. 10, filled, color = "Yellow"):

B := plots[animate](plot, [[2*r*cos(phi)/Pi, 2*arcsin(sin(r*sin(phi)))/Pi, phi = 0 .. 2*Pi], thickness = 3, color = blue, scaling = constrained], r = 0 .. 5*Pi, frames = 100):

plots[display](A, B);

 Addition - animation in 3D:

A := plot3d([x,cos(t),sin(t)], x = -10 .. 10, t=0..2*Pi, style=surface, color = "Yellow", transparency=0.75):

B := plots[animate](plot3d, [[2*r*cos(phi)/Pi, 2*arcsin(sin(r*sin(phi)))/Pi*cos(t), 2*arcsin(sin(r*sin(phi)))/Pi*sin(t)],phi = 0 .. 2*Pi, t=0..2*Pi, style=surface,color = "Niagara Azure", scaling=constrained, numpoints=10000], r = 0 .. 5*Pi, frames = 100):

plots[display](A, B, axes=none);

Better quality of plotting without additional option  grid  enables the use of spherical coordinates.

Compare:

restart;

A:=plots[implicitplot3d](x^2+y^2+z^2-2^2, x = -2 .. 2, y = -2 .. 2, z = -2 .. 2, style = surface, axes = normal, view = [(-2.7 .. 2.7)$3],  grid = [30, 30, 30]):

B:=plot3d([2*cos(s)*sin(t), 2*sin(s)*sin(t), 2*cos(t)], s=0..2*Pi, t=0..Pi, style=surface):

plots[display](<A | B>, axes = normal, view = [(-2.7 .. 2.7)$3]);

Use multiple assignments as in the example below:

a, b := solve(x^2-1=0);

a;

b;

Another way:

plots:-display( plot3d(exp(x*y), x= 0..1, y= 0..1),   plots:-intersectplot(y=exp(x), z=0, x= 0..1,

y=0..exp(1), z=0..exp(1), color= red, thickness= 2, axes=normal));

restart;

p := proc()

local A, B, D, m, n, k, r;

A := Matrix([[0, 1, 0, 0, 0], [1, 0, 1, 0, 0], [0, 1, 0, 1, 0], [0, 0, 1, 0, 1], [0, 0, 0, 1, 0]]);

  r := LinearAlgebra:-RowDimension(A);

  D := Matrix(r); 

   for k from 1 to r-1 do

     B := A^k; 

       for m from 1 to r do

           for n from 1 to r do

               if m <> n and B[m,n] <> 0 and D[m,n] = 0 then D[m, n] := k;

               end if;

           end do;

        end do;

   end do;

   D;

end proc:

p();

Tools->Options->Interface->Default zoom