Follow these steps:

1. Sketch a graph of the function defined in your equation (6).
2. Did you do that?
   If not, then go to Step 1.
3. You DID sketch a graph of the function defined in equation (6), right?
   Now go to step 4.
4. What is the area under that graph?  You should be able to tell that the area is zero by just looking at it.  No computation is necessary.
5. Done!

Aside: Where you say Kronecker delta you should say Dirac delta.  Kronecker delta is something quite different.  What you have is neither.

In your worksheet you have:

x := seq(x =-1..1, dx);

That doesn't make much sense.  Change that to:

x := Array(0..nx, i-> 2/nx*i - 1, datatype=float[8]):

Afterward, when you want to refer to the entries of x, use square brackets, as in x[i].

There is no solution in terms of elementary functions, but you can obtain a solution numerically if you provide intitial conditions, as in:

de := -1.699678499+diff(z(t),t,t) - 3.813322790*z(t) + 0.6966019659*z(t)^2 + 0.01548004368/z(t)^2 = 0;
ic := z(0)=1, D(z)(0)=0;
dsol := dsolve({de,ic}, numeric);
plots:-odeplot(dsol, t=0..8);

which produces the following plot:

Note that in specifying your differential equation, every z should be entered as z(t).

Have a look at this worksheet.  It does some calculations with your problem but I don't quite understand what you are looking for, so it may not be quite finished yet.

partial-solution.mw

Maple is attempting to produce an answer to an ill-posed problem.  Neither answer satisfies your initial condition D(y)(0) = 0.  You can see that by plotting them:

plot((1/3)*sin(3*t), t=0..2);
plot((1/3*(Heaviside(t)-1/2))*sin(3*t), t=0..2);

That's not a fault of Maple's. Specifying a delta function acting at t=0 is not meaningful. Consider applying it at a later time, as in:

eq := diff(y(t),t,t) + 9*y(t) = Dirac(t-1);

Remove the assignments x:=... and  y:=....  Instead, plot the result through:

odeplot(sol, [u(s)*sin(varphi(s)), -u(s)*cos(varphi(s))], s=0..10);

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.