Ronan gave a link to my post from almost two years ago where I illustrated the use of quaternions in tracking rotations in space.  That motivated me to clean up the code a bit and then apply it to the phenomenon of the flipping T-handle which is the subject of the current thread.  The code, which combines Euler's equations of motion and quaternions, is sufficiently commented to let you adapt it to other situations.

Here is the worksheet T-handle-flip.mw and here is the resulting animation:

@ecterrab You are correct in stating that this problem has "no solution".  To be precise, this problem has no classical solution.

A great deal of effort was devoted throughout the 20th century to make sense of exactly such boundary value problems which arise in just about every application, especially in physics.  The concept of generalized derivatives and weak solutions were introduced precisely for that purpose.  If we dismiss weak solutions, we will be throwing away most of the useful applications of the PDEs.

Consider, for instance, a slender (one-dimensional) bar currently at a uniform temperature of 100 degrees.  We bring one end of it in contact with an ice cube.  We want to determine how the bar's temperature evolves in subsequent times.

The temperature, u(x,t), satisfies the PDE u_t = u_xx, where subscripts denote differentiation.  The initial condition is u(x,0)=100.  The boundary condition is u(0,t)=0.  Certainly these two conditions are incompatible but we don't want throw up our hands and say that this problem is not solvable.  True, no classical solution exists, but a weak solution does, and that's what's relevant to this question.

Mathematica's solution that was posted earlier (and is now deleted) expresses the solution of the PDE in the Fourier series.  The Fourier series solution picks up the weak solution, therefore that solution was correct in the sense that it produced the expected weak solution to the problem.  In the worksheet below I show how to reproduce Mathematica's solution with Maple's pdsolve().  The process hits a few snags but eventually we arrive at the desired solution.

P.S.: To be clear, I am not here to berate the terrific job that you have done (and continue doing) with extending Maple's capabilities in many directions.  All I want to point out is that the Fourier series picks up the weak solution, and that it should not be dismissed.

[Worksheet produced in Maple 2018.1]

Download pde-bc.mw

I don't quite understand what you mean by "link a worksheet written in Maple 16 into WordPress" but let me try.

WordPress can display HTML code as well as images.  In Malpe's menubar select File->Print and save your worksheet as a graphics image.  On my Maple 2018 on Linux the image is saved in the PostScript format.  In your Maple 16 it may be something else, but it does not matter much.  You may convert one image into any other image format such as PNG, which then you pass on to WordPress.

plot([0, sin(x)], x=0..2*Pi, axes=boxed, color=[black,blue]);

In the Layout palette pick the entry that consists of a green "A" on top of a red "b".  Then edit to replace the "A" with "max" and the "B" with the desired subscript.  Here is what we get:

This works in Maple 2018.  I don't have Maple 2015 to test.

Download worksheet: mw.mw

The answer depends on what you assume about x.

Is x itself bounded?  If yes, then cos(i*x) is bounded.

Is x imaginary and unbounded?  If yes, then i*x is real, and therefore cos(i*x) is bounded because sin(i*x)^2 + cos(i*x)^2 = 1.

Is x real and unbounded? If so, then:

convert(cos(I*x), exp) assuming x::real;

shows that cos(i*x) is 1/2 exp(x) + 1/2 exp(-x), and therefore cos(i*x) is unbounded.

And here is an exercise for you: Suppose x is complex.  What can you tell about the boundedness of cos((i*x)?

The purpose of phaseportrait is to plot solutions of a system of differential equations.  There is no differential equation in what you have shown.

Suggestion: Maple's help page on phaseportrait is quite limited and not very instructive.  To its credit, it refer the reader to DEplot's help page which does what phaseportrait is supposed to do, and is much better documented.  Look up the examples in DEplot's help page and come back here if you still have questions.

The abs(t) seems to confuse Maple. The following equivalent alternative works:

Download mw.mw

Download mw2.mw

Let me write L instead of LL for greater legibility.  By differentiating your expression for L we get:

Solve this as a differential equation:
dsol := dsolve({de, ic}, numeric);

Then plot it:

plots:-odeplot(dsol), phi=0..2*Pi);

To read of the value of the solution at phi=4, do:

dsol(4);

iR := rand(1..nR):
iC := rand(1..nC):
C(iR(),iC());  # pick a random element of C
C(iR(),iC());  # pick another random element of C

To define the matrices A[k]:

assign(seq(A[k] = <k,k^2;2*k,2*k^2>, k=1..10));

To make a block-diagonal matrix out of these:

M := LinearAlgebra:-DiagonalMatrix([seq(A[k], k=1..10)]);

To view M:

evalm(M);

Aside: If you have no need to refer to the individual matrices A[k], then you may bypass the middleman and make the block-diagonal matrix directly:

M := LinearAlgebra:-DiagonalMatrix([seq(<k,k^2;2*k,2*k^2>, k=1..10)]);

The differential equation itself expresses y''(t) in terms of y(t) and y'(t).  Therefore once you have obtained y(t) and y'(t), you can evaluate y''(t) directly, like this.

Download mw.mw

Regarding your statement:
I know function f(x) for y>1.4266 is monotically decreasing and there is no turning points,(actually i can compute infection points of above function: y=1.4266 and x=4/5 ).

That's correct but the calculation in your worksheet may be simplified to one line:
solve({diff(F(x,y),x) = 0, diff(F(x,y),x,x) = 0});

As to your statement:
I know function for 1.4266>y>1.0144 has two turning points a minimum for x<4/5 and a maximum for x>4/5.

I don't see the significance of the 1.0144.  The function has a minimum and a maximum for all values of y in 0 < y < 1.4266.  You can see that by letting
F := (x,y) -> x^2+y*(1-x)^(3/4);
and then plotting F(x,y) for 0 < x < 1 and various values of y.