Convert the system to a system of first order equations:

restart; with(DEtools):
M:=1: L:=1: g:=9.81:
sys := diff(x(t),t) = dx(t), diff(z(t),t) = dz(t),
diff(dz(t),t)=-1/z(t)*diff(x(t),t)^2-1/z(t)*diff(z(t),t)^2-x(t)/z(t)*diff(dx(t),t),
M*(-L^2/(x(t)*z(t)))*diff(dx(t),t)=M*g-M*(-1/z(t)*diff(x(t),t)^2-1/z(t)*diff(z(t),t)^2);

Now use DEplot to plot z vs x:

DEplot( {sys}, [x,z,dx,dz], t=0..2,
     {[dx(0)=0,x(0)=0.99999999,dz(0)=0,z(0)=sqrt(L^2-0.99999999^2)]},
     stepsize=0.001, scene=[x,z], x=-2..2, z=-.1..1.1, linecolor=blue );

DEplot.mw

Of course you can calculate the eigenvalues, and check if thesse are all positive. In this case these are rather complicated expressions. Easier - especially when your matrices become larger - is to calculate the determinant of submatrices. If these are all positive, the quadratic form is positive definite; if the signs of these determinants alternate, it is negative definite. So, for your example:

restart;
A:=Matrix([[8,4,2,1],[4,8,2,1],[2,2,8,1],[1,1,1,8]]):
seq( LinearAlgebra:-Determinant(A[1..i,1..i]), i=1..4 );
                        8, 48, 352, 2736

All positive, so the quadratic form is positive definite.

Your system is underdeterminated, zo there are many free parameters. A naive approach is to set all free parameters zero. I give a small example:

restart; with(LinearAlgebra):
A := RandomMatrix(4,9);
b := <0,3,0,0>;
c := LinearSolve(A,b);
eval(c,indets(c)=~0);

Another possibility is:

dsnx := dsolve(sysx, numeric, method = bvp[midrich], continuation = lambda, 
       maxmesh = 20000, initmesh = 20000, abserr = 10^(-2), 
       output = Array([seq(i,i=1..10,0.1)]) );
res := dsnx[2,1]:
res[[1,10,20,40,50,90],[1..3]];

I added de derivative of the third equation, and I make all variables functions of t:

restart;
a := miu1(t) = Diff(lambda1(t),t) + lambda3(t);
b := miu2(t) = Diff(lambda2(t),t) + lambda1(t) - 4*lambda3(t);
c := miu3(t) = Diff(lambda3(t),t) + lambda2(t) - 3*lambda3(t);
c1 := Diff(lambda3(t), t, t) = -(Diff(lambda2(t), t)) + 3*(Diff(lambda3(t), t)) + Diff(miu3(t), t);
d := miu(t) = 2*Diff(lambda3(t),t$2) - 5*Diff(lambda3(t),t) + 7*lambda3(t);

Now I use the first three equations to solve for the derivatives of lambda, and substitute rhe result in the last equation:

 s := solve( {a,b,c}, {Diff(lambda1(t),t),Diff(lambda2(t),t),Diff(lambda3(t),t)} ):
A1 := subs(s,d);

This contains the second derivative of lambda3, which is also in the third equation. But this equation, in turn, contains again derivatives of lambda. So I use the 2nd ant 3rd equation to get expressions for these derivatives in terms of lambda:

A2 := solve( {b,c}, {Diff(lambda2(t),t),Diff(lambda3(t),t)} );
A3 := subs(A2,c1);
subs(A3,A1);

I don't understand your second question.

Maple's answer is correct:

restart;
solve([-sin(x) = cos(x), 2*Pi > x and x > 0], x, AllSolutions = true);
                     /      1            \ 
                    { x = - - Pi + Pi _Z1 }
                     \      4            / 

about(_Z1);
Originally _Z1, renamed _Z1~:
  is assumed to be: AndProp(integer,RealRange(1,2))

This shows that _Z1=1 or _Z1=2.

DE := diff(y(x),x$2) - 3*diff(y(x),x) - 4 = exp(x)*sin(4*x);
dsolve(DE);

Try it yourself the next time

normalcdf:= proc(x1,x2,mu,sigma)
local t, dist:
  dist:= t -> evalf(Statistics:-CDF(Normal(mu, sigma), t));
  print(P(z<x1)=dist(x1));
  print(P(z>x2)=1-dist(x2));
  print(P(abs(z-mu)>sigma)=1-dist(x2)+dist(x1));
  print(P(abs(z-mu)<sigma)=dist(x2)-dist(x1))
end proc:

normalcdf(2100,2200,2150,50);
                   P(z < 2100) = 0.1586552540
                   P(2200 < z) = 0.1586552540
               P(50 < |z - 2150|) = 0.3173105080
               P(|z - 2150| < 50) = 0.6826894920

Notice that your function has two zeros at x=0 for all omega, so divide by x^2. To plot the other zero, depending on omega:

restart;
f := x + omega^2*x^3 - omega:
pnts := [seq( [omega, fsolve(f)], omega=0..10,0.1 )]:
plot( pnts );

I suppose that you can write a Newton-Raphson procedure to use instead of fsolve if it is essential for your problem to use this method.

Your second question is to check the solution:

restart:
g:= x-> d*x^2+e*x+f:
sol := solve({g(0)=4, g(3)=7, g(-2)=18}, {d,e,f}):
G := subs( sol, eval(g) );
     8  2   19      
x -> - x  - -- x + 4
     5      5       
is( G(0)=4 and G(3)=7 and G(-2)=18 );
                              true

Replace 2sin(pix)  by 2*sin(Pi*x) : Pi with capital P, because pi is only a name for a variable.

If the list of x-values is called X, and the lists of y-values are called Y1,...Ym:

L := [seq( seq([X[j],Y||j[i]], i=1..n), j=1..m )]:
plots:-pointplot(L);

replace lengte := 22 by lengte := 22.0

Your notation for f suggests that you mean a 3d vectorfield in cartesian coordinates. See ?VectorCalculus

Using this package, the unit vectors i, j, k are represented as e_x, e_y, e_z.

You will get a nicer plot if you make a parametrization of the surface:

plot3d( [13.2*sin(phi)*cos(theta), 24.3*sin(phi)*sin(theta), 22*cos(phi)], 
             phi=0..Pi, theta=0..2*Pi, scaling=constrained );