At the end of your worksheet insert:

series(rhs((3)), b);

For illustration I have made up somewhat arbitrary functions f, g, h, and k, and consequently the plots are not particularly pretty, but the method is general; replace f, g, h, and k by functions of your choice.

restart;
with(plots): with(plottools):
k := (x,y,z) -> x*sin(y) + y*sin(z) + z*sin(x); 
p := implicitplot3d(k(x,y,z)=0, x=Pi..3*Pi, y=Pi..3*Pi, z=-Pi..Pi, grid=[30,30,30]):
display(p, orientation=[-25, 75, 0 ]);

f := (x,y,z) -> x;
g := (x,y,z) -> y;
h := (x,y,z) -> sin(z);
T := (x,y,z) -> [ f(x,y,z), g(x,y,z), h(x,y,z) ];

Your request amounts to asking what the mapping T does to the surface p.  Here it is:

transform(T)(p);

Or, after applying some cosmetic touches:

display(transform(T)(p), style=surface, color=gold,
    orientation=[15,70,0], scaling=constrained);

The standard arcsin function expresses the inverse of the sin function over the range [-Pi/2,Pi/2].  In your function, the arcsin function needs to be interpreted more broadly.  To remove the ambiguity, it is best to get rid of the arcsin.

Thus, in the equation
     arcsin(sqrt(S))-sqrt(S)*sqrt(1-S) = t
isolate the arcsin(sqrt(S)) term, then apply sin to both sides of the equation.  Here is how:

arcsin(sqrt(S))-sqrt(S)*sqrt(1-S) = t;
isolate(%, arcsin(sqrt(S)));
myeq := sin(lhs(%)) = sin(rhs(%));

Now plot:

plots:-implicitplot(myeq, t=0..Pi, S=0..3,
     grid=[20,20], gridrefine=2);

Remark:
I must add that implicitplot should be a solution of last resort.  Once you see what is going on, you may produce the same plot without implicitplot as follows:

Arcsin := S -> Pi - arcsin(S);
U := arcsin(sqrt(S))-sqrt(S)*sqrt(1-S);
V := Arcsin(sqrt(S))-sqrt(S)*sqrt(1-S);
plots:-display([
  plot([U, S, S=0..1]),
  plot([V, S, S=0..1])
], labels=[t,S]);

The teaching method that you are proposing would have been appropriate if the students did not have Maple.  In the presence of Maple, your students may find that approach puzzling.

What is the point of calculating the slope to the left and to the right of an extreme point when Maple can easily calculate the slope everywhere. A graph of f(x) and the derivative f'(x) can express your idea much more forcefully.

f := x -> x^3-x^2-4*x+12;
g := D(f);
solve(g(x)=0);
plot([f(x),g(x)], x=-2..3, color=[red,blue]);

Here is something to get you started:

mw.mw

Almost any textbook on elementary differential equations has a chapter on the Laplace transform.  That provides the simplest way to solve your system of equations although the final form of the solution itself is quite messy.

If you don't know about the Laplace transform, you should pick up a book and read about it.  You won't learn about it by reading this post.

If you do know about the Laplace transform, then you should try applying it to the simple equation

$\mathrm{de}≔\frac{d^2}{d{t}^2}x\left(t\right)+bx\left(t\right)=f\left(t\right)$

The file mw1.mw shows how to do that in Maple.

Once you understand how that works, then you should examine the file mw2.mw which shows how to solve your system of equations.

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]);