@Markiyan Hirnyk

1. The inequality  (1/2)*(x-y)^2-(x-y)^4>=0  follows from the equation  (1/2)*(x-y)^2-(x-y)^4 = (-2*x^2+y^2)^2

2.  I also was solving it as  @_Maxim_   but in Maple 2016  I got this:

solve({y >= 4*x^4+4*x^2*y+1/2, (1/2)*(x-y)^2-(x-y)^4 = (-2*x^2+y^2)^2, -2*x^2+y^2>=0}, [x, y]);

                                                           [ ]

3. In any case, visualization is useful, since it allows you to test the correctness of the solution visually. You do not require Maple any proof! How do you know that the solution received is correct?

@Markiyan Hirnyk Yes you are right, I incorrectly determined those regions in which the solution is located. Here is the right picture:

A:=plots:-implicitplot([y = 4*x^4+4*x^2*y+1/2, (1/2)*(x-y)^2-(x-y)^4 = (-2*x^2+y^2)^2, -2*x^2+y^2=0], x=-3..3, y=-3..1, color=[red,blue, green], thickness=2, gridrefine=5):
B:=plots:-inequal({-2*x^2+y^2>=0, y>=4*x^4+4*x^2*y+1/2}, x=-3..3, y=-3..1, color="LightGreen"):
plots:-display(A, B, scaling=constrained);  

So the final result is 2 points:  [x = -1, y = -3/2], [x = 0, y = 1/2]

@Jason Lee  You wrote "...can u explain to me how to use implicitplot in Maple." Preben already showed how to use this command. I can only add that the roots of the system of equations are only the intersection points of curves corresponding to different equations of the system. Therefore, in the plot, you should use different colors for different curves:

eqn1:=-x*y^2+4*x=5;
eqn2:=(1/3)*x^3+y^2=1;

plots:-implicitplot([eqn1, eqn2], x=-10..10, y=-10..10, color=[red,blue], gridrefine=3);

@MrYouMath  You did not copy, but just reprinted the code, because should be  symbol  instead of  symobl .

@MrYouMath  No, Maple does not understand such expressions as  {p__11,p__12,...p__1n, p__12,...,p__2n,...,p__1n,....,p__nn}. You must write this code in the form of a procedure

restart;
LyapunovEq:=proc(A::Matrix)
local n, P, Id, eqn;
uses LinearAlgebra:
n:=op([1,1], A);
P := Matrix(n, symbol=p);
Id := Matrix(n, shape=identity);
eqn := Transpose(A).P+P.A =~ -Id;
solve({seq(seq(eqn[i,j], j=1..n), i=1..n)}, indets(P));
end proc: 

Examples of use:

LyapunovEq(<1,2; 3,4>);
LyapunovEq(Matrix(3, symbol=a));

 

@Carl Love Yes, but a little later. I already wrote that this problem is only part of an entertaining problem. I prepare its complete solution  as a separate post, it remains to finish some details. My algorithm for selecting the corresponding matrices will simply be part of this post (not to repeat it twice).

@Carl Love  Thank you!

@Carl Love Thank you very much. Impressive method and its improvement in the comment "Under 1 second". Unfortunately, it seems to work only if each row and each column has exactly 2 zero or  2 units as in the original problem.

@vv Thank you. Unfortunately your code does not work on my computer for some reason:

Only with the option  compile=false  I get the result, but the time is 10 times longer.

When I wrote about my method for 100 seconds, I used another technique (uncompiled code). I'll submit this code later.

@tsunamiBTP  I did this and did not find any fundamental discrepancies. Some differences are due to the fact that in  M_new  I am greatly rounding the results (evalf[3]). For comparison, I calculated the hundredth line of my and your matrix. The graphs are indistinguishable. See untitled5_new2.mw

@John Fredsted  About this construction see Matrix and Vector Construction Shortcuts help page. By the same method, we can remove from a matrix any row or several rows or columns:

A:=LinearAlgebra:-RandomMatrix(6, generator=0..9);
A1:=<A[..3], A[5..]>;  # 4th row deleted of A
A2:=<A[..,..4] | A[..,6]>;  # 5th column deleted of A

I think this method is faster than calling  LinearAlgebra:-DeleteRow  or  LinearAlgebra:-DeleteColumn .

@tsunamiBTP  See corrected file  untitled5_new1.mw 

@tsunamiBTP  You wrote "This seems to work like a nested FOR loop". No, these are the different things. Here is a simple example of use nested seq. Let's say you want to create a matrix 10x10 that has symbols  x1, x2, ... ,x10  in the first line, the second line contains the same symbols, which squared, then cubed and so on:

Matrix([seq([seq((x||j)^i, j=1..10)], i=1..10)]);

Another example. It is very convenient and effective to use  seq  command  to obtain the sequence or list of objects with specified properties. In the example, we generate the sequence of Pythagorean triplets with  hypotenuses <=100:

seq(seq(seq(`if`(a^2+b^2=c^2,[a,b,c],NULL), b=a+1..c-1), a=1..c-2), c=3..100);

[3, 4, 5], [6, 8, 10], [5, 12, 13], [9, 12, 15], [8, 15, 17], [12, 16, 20], [7, 24, 25], [15, 20, 25], [10, 24, 26], [20, 21, 29], [18, 24, 30], [16, 30, 34], [21, 28, 35], [12, 35, 37], [15, 36, 39], [24, 32, 40], [9, 40, 41], [27, 36, 45], [14, 48, 50], [30, 40, 50], [24, 45, 51], [20, 48, 52], [28, 45, 53], [33, 44, 55], [40, 42, 58], [36, 48, 60], [11, 60, 61], [16, 63, 65], [25, 60, 65], [33, 56, 65], [39, 52, 65], [32, 60, 68], [42, 56, 70], [48, 55, 73], [24, 70, 74], [21, 72, 75], [45, 60, 75], [30, 72, 78], [48, 64, 80], [18, 80, 82], [13, 84, 85], [36, 77, 85], [40, 75, 85], [51, 68, 85], [60, 63, 87], [39, 80, 89], [54, 72, 90], [35, 84, 91], [57, 76, 95], [65, 72, 97], [28, 96, 100], [60, 80, 100]

@Markiyan Hirnyk  Here is a check (I just added a few commands to your document):

Download Check.mw

@mehdibaghaee Sorry, I did not notice the second condition.
Should be:

seq(`if`(V[i]>0 and V[i]<>infinity, [i, V[i]], NULL), i=1..10);