The boundary condition D[1](H,0)=0 does not make sense.  Probably you meant D[2](H,0)=0.

But after fixing that, you will get another error message:

     Error, (in pdsolve/numeric) unable to handle elliptic PDEs

That's because Maple's PDE solver can handle only limited types of PDEs.  Your equation is not among them.  You may consider writing your own solver (implementing a finite difference approximation scheme is not too difficult) or finding another software that can do this for you.

The purpose of a differential equation is to define the derivative y' in terms of x and y.  Your differential equation is meaningless at x=1 because it reduces to y' - 2y'^3 = 0 which has three real roots:
   y' = 0,   y' = 1/sqrt(2),    y' = -1/sqrt(2).
Consequently, at x=1 your differential equation is unable to tell you which way to go.

The x=1 is not very special.  That same issue exists at any x where the algebraic equation y' - 2y'^3 = -2x + 2 has multiple real solutions y'.

You may easily verify that the equation y' - 2y'^3 = -2x + 2 has a single real solution y' if x is in the interval x=0 to x=b, where b = 1-sqrt(6)/18 = 0.8639...

On that interval you may solve and plot the solution of your initial value problem.

Details are in the attached worksheet.

mw.mw

I think by "phase transition" you mean "phase plot".

Look up Maple's help on DEplot3d.  There are several instructive examples there.

The VariationalCalculus package is made for this purpose.

In the example you have shown, I assume that l is l(t) and ld is its derivative.  Therefore this is how we proceed.

restart;
with(VariationalCalculus):
T := 1/6*rho*diff(l(t),t)^2*(2*L^2 + 4*H^2/l(t)^6);

ans := EulerLagrange(T, t, l(t));

The "ans" above consists of two parts.  The first part is your problem's Euler-Lagrange expression.  Set it equal to zero to get the Euler-Lagrange equation.  The second part is the integrated version of the Euler-Lagrange equation, known as the first integral.  You may extract one or the other part as follows.

Extract the Euler-Lagrange expression:
remove(type, ans, 'equation');

Extract the first integral:
select(type, ans, 'equation');

restart;
plot([2*sqrt(cos(2*t)), -2*sqrt(cos(2*t))], t=-Pi/4..Pi/4, coords=polar, color=red);

Your calculations need a couple of fixes.

1. The subs() command doesn't do what you mean.  To see the problem, do
   subs({G=1,M=1},satODE1);    # BAD
   To fix, change that to:
   subs({G=1,M=1},[satODE1]);  # GOOD
2. Your ics asks Maple to plot four solutions, but that's not what you want.  You want one solution with initial conditions specified for the four variables.  To fix, change ics to [ics] in your plotting command.

Applying these fixes we get:

mysys := subs({G=1,M=1},[satODE1]);

DEplot(mysys,{x(t),y(t),vx(t),vy(t)},t=-2..2,[ics],scene=[x(t),y(t)],scaling=constrained);

I assume you use Maple's palettes to produce the accents.  An accented character such as $\overrightarrow{\beta }$ is produced through the code

`#mover(mi("β"),mo("→"))`

where → is the standard HTML code for "right arrow".

Then it would have made sense if Maple's palettes produced $\stackrel{—}{\beta }$ through

`#mover(mi("β"),mo("—"))`

but unfortunately  it doesn't.  Instead, it pruduces $\stackrel{—}{\beta }$ through conjugate(beta), where conjugate is Maple's standard operator for the complex cojugate.  That works fine for undefined symbols but not for symbols with numeric values.  For instance you cannot make a "5" with a bar over it because the complex cojugate of a real number is itself.  I consider the use of conjugate for this purpose a bad design decision, if not a bug.

Here is where the trouble with gamma comes in.  The symbol gamma is predefined in Maple to Euler's constant.  Do evalf(gamma) to see it.  As a result, conjugate(gamma) is just gamma, that's why it is not shown with a bar.

So, how do we make gamma with a bar?  Here is how.  Let

macro(bgamma = `#mover(mi("gamma"),mo("__"))`);

Then bgamma will produce $\stackrel{—}{\gamma }$.

Probably you are thinking that A, B, alfa, and y are real, but Maple doesn't know that.  You have to tell it what's in your mind.

Additionally, the "i" in the second term should be "I".

Additionally, probably by "alfa" you mean "alpha".

Fixing these, we get:

w := A*exp(-alpha*(1+I)*y)+B*exp((1+I)*y);

u := Re(w) assuming A::real, B::real, alpha::real, y::real;

See if this meets your needs.  I am plotting simple functions here ( 0.5*exp(-x) and exp(-x) ) for illustration.  Replace them with yours.

restart;
with(plots): with(plottools):
p1 := plot([0.5*exp(-x), exp(-x)], x=0..1, y=0..1, color=[blue,"Green"]);
T := (x,y) -> piecewise(y>x, [(x+y)/2,(x+y)/2], [x,y]):
p2 := transform(T)(p1);
p3 := plot(x, x=0..1, color=red, thickness=2);
display([p2,p3], view=[0..1,0..1]);

I see that you have been struggling with animation of rigid body motion.  I think a good deal of the trouble is due your dependence on the ugly/messy Euler angles and Euler rotation matrices.  If you move forward from Euler's time by about a century, you will discover the beautiful algebra of quaternions and their use in describing rotations in 3D.

Unfortunately these have not found their way yet into the standard curriculum.  At one point I started writing a textbook on the subject based on my lecture notes, but have not had the time/energy/motivation to finish.  Maybe some day…

In the mean time, you may find the attached Maple worksheet of some use.  It illustrates the use of quaternions in rigid body dynamics.  You will see that there are no Euler angles or rotation matrices in it.  I think the contents should be easy to grasp and modify to suit your needs although it gives no explanation of the theory behind it.  Here is a sample animation produced in that worksheet:

rigid-body-dynamics-with-quaternions.mw

Edit:
But wait!  Looking over this now I see that in my attemp to simplify things I may have gone too far and written nonsense.  Don't bother with it for now.  I will double-check and upload a corrected version later.

Edit:
OK, false alarm. What I had posted was correct after all. Here I am posting a revised version which expands some of the documentation.  The calculations remain the same.

rigid-body-dynamics-with-quaternions-v2.mw

Note that the pair of algebraic equations 2A + B = 0, 2B + A = 0 implies that A = 0 and B = 0.

Comparing these with your eq2 and eq3, we conclude that the second derivatives of y(x) and z(x) are zero, and therefore their first derivatives are constants.

But that contradicts your boundary conditions D(y)(0)=0, D(y)(1)=1, hence no solution.

The boundary condition D[1](u)(0,t)=1000 should be D[1](u)(0,t)=-1000.  Remember that heat flows from hot to cold, that is, against the gradient.

Method 1:

pointplot([seq([n/L, g[n]], n = -L .. L)],  color = blue, connect);

Method 2:

plot([seq([n/L, g[n]], n = -L .. L)],  color = blue);

The equation is solvable if and only if

To see this, forget about the boundary conditions for the moment.  Find (by hand, not Maple) the general solution of the differential equation which will involve two arbitrary constants, let's say C₁ and C₂.   This is done through what is known as the method of variation of constants which you will find in any elementary differential equations textbook.

Then applying the boundary conditions leads to two linear algebraic equations in the two unknowns C₁ and C₂.  You will immediately see that the solvability of the system is equivalent to the integral condition stated above.

Advice: It is best to do this in paper and pencil first if you want to understand what is going on. Afterward you can do the same thing in Maple (it's very easy).  If you get stuck, show your work then I or someone else will help you.

<(f(t) $t=-3..3)>;