To find all the roots of an analytic function in a certain range it is better to use the  RootFinding[Analytic]  command:

v := unapply(H*sin(w*t), t):

L := 0.080:

H := 0.020:

Vf := 0.3:

w := 10:

zeros1 := sort([RootFinding[Analytic](v(t), t, re = 0.5 .. 2, im = -1 .. 1)]);

zeros2 := sort([RootFinding[Analytic]((D(v))(t), t, re = 0.5 .. 2, im = -1 .. 1)]);

plot([seq([zeros1[i], t, t = -0.01 .. 0.01], i = 1 .. 5), v(t)], t = 0.5 .. 2, color = [yellow$5, blue], linestyle=[3$5,1], thickness=[2$5,1]);

Similarly, you can build  vertical lines at the points of extrema.

in Maple 2015.1

z := x+I*y:

Expr := evalc(abs((1/2)*z^2+3*z-1));

plots[inequal](Expr <= 1, x = 0 .. 0.8, y = -0.4 .. 0.4, scaling = constrained);

Addition: It is better to paint the inside of the region and its border not by the default but using the appropriate options explicitly. Since in this example, we have an open region, its boundary is shown by a dashed line: 

z := x+I*y:
Expr := evalc(abs((1/2)*z^2+3*z-1)):
plots[inequal](Expr <= 1, x = 0 .. .8, y = -.4 .. .4, optionsfeasible = [color = yellow], optionsclosed = [color = blue, linestyle = 6, thickness = 2], scaling = constrained);

The added vectors does not necessarily lie in the orthogonal complement to the linear manifold of  initial rows. Procedure  P  tests given matrix  A , checking the conditions " given a matrix M[n*2n] of n linearly independent row vectors" and adds to it the missing  n  linearly independent unit vectors.

P:=proc(A::Matrix)

local n, m, A1;

uses LinearAlgebra;

n:=RowDimension(A); m:=ColumnDimension(A);

if m<>2*n then error "The column dimension should be 2 times greater than row dimension" fi;

if Rank(A)<n then error "The rows of matrix are linearly dependent" fi;

A1:=<A; IdentityMatrix(m)>;

Basis([seq(A1[i],i=1..n+m)]);

<op(%)>;

end proc:

Your example:

A:=<4,1,0,0; 0,0,4,1>;

P(A);  

LinearAlgebra[Rank](A); 

A:=<4,1,0,0; 0,0,4,1>;

V1:=<0|0|0|1>; V2:=<0|0|1|0>;

<A; V1; V2>;

I don't know what is rref form. 

I displayed the results of the last few commands. It seems determinant  ZZ  calculated correctly. But obviously it has no real roots, so you get an empty list  Sol:=[ ] . Therefore the following command gives infinity:

ZZ := Determinant(q);

Sol := [fsolve(ZZ, omega)];

J := min(select(`>`, Sol, 0));

OMEGA := min(select(`>`, Sol, 0));

P2 := subs(omega = OMEGA, Q2);

Separately:

min(select(`>`, [ ], 0));

                           infinity

Just specify both surfaces by implicit equations (in Cartesian coordinates), and then  plots[intersectplot] :

plots[intersectplot](x^2+y^2+z^2 = 1, sqrt(3)*x-y = 0, x = 0 .. 1, y = 0 .. 1, z = -1 .. 1, color = red, thickness = 2, scaling = constrained, axes = normal);

 You can also directly build in cylindrical coordinates:

plots[intersectplot](r^2+z^2 = 1, phi = (1/3)*Pi, r = 0 .. 1, phi = 0 .. (1/2)*Pi, z = -1 .. 1, color = red, thickness = 2, scaling = constrained, axes = normal, coords = cylindrical);

Edited.

After  assign(%):  you can not write  x = 0.1e-2 .. 10 . Should be  'x' = 0.1e-2 .. 10

Corrected version  (in a classic worksheet Maple 2015.1):

restart;

with(plots):

eq1 := z = y*log(x): eq2 := z = y+x*log(x): 

DispIntersecting := implicitplot3d([eq1, eq2], x = 0 .. 10, y = -30 .. 30, z = -40 .. 40, color = [blue, green]):

solve({eq1, eq2}, [x, y, z]):

assign(%):

DispIntersection := spacecurve([x, y, z], 'x' = 0.1e-2 .. 10, color = red, thickness=3, view = [0 .. 10, -30 .. 30, -40 .. 40]):

display(DispIntersecting, DispIntersection, axes = boxed, scaling = constrained, labels=['x','y','z']);

restart;

for k to 10 do

Sol := dsolve({diff(x(t), t, t)+k*sin(x(t)) = 0, x(0) = (1/6)*Pi, (D(x))(0) = 0}, numeric):

f := t->rhs(Sol(t)[2]):

A[k] := seq(plottools[line]([0, 1], [sin(f(t))/k, 1-cos(f(t))/k], thickness = 2), t = 0 .. 24, 0.2):

B[k] := seq(plottools[disk]([sin(f(t))/k, 1-cos(f(t))/k], 0.012-(k-1)*0.0006, color = black), t = 0 .. 24, 0.2):

F[k] := seq(plots[display](op([A[k][i], B[k][i]])), i = 1 .. 121):

od:

plots[display](seq(plots[display](op([seq(F[k][i], k = 1 .. 10)])), i = 1 .. 121), insequence = true, scaling = constrained, axes = none);

You do not need any temporary variables. You can just do this by one command

subs([b=a, c=b, a=c], a*b*c+a*b);

                       a*b*c+a*c

See help on  subs  for details.

Corrected version:

dsol1 := dsolve({diff(f(eta), eta, eta, eta)+3*f(eta)*(diff(f(eta), eta, eta))-2*(diff(f(eta), eta))^2+1 = 0, f(0) = 0, (D(f))(0) = 0, ((D@@2)(f))(1) = 0}, numeric, output = procedurelist);

plots[odeplot](dsol1, [eta, f(eta)], eta = 0 .. 1, color = red);

Addition:  in your  2D-math replace  2  in the boundary condition for the second derivative by  (2) 

Maybe OP meant the creation of new functions from existing functions by using the composition of the operators and other operators?

An example:

x:=t->t^2:   y:=t->sin(t):

f:=x@(y+x):

g:=y@(x-y):

f(t);

g(t);

To find the values of a function better  to specify it as a function not as an expression:

g:=y->200*y:

g(4);

              800

The workaround for your method:

g(0.5,y):=200*y:

eval(g(0.5,y), y=4);

               800

From the standpoint of rigorous mathematics, evalf  command is not sufficiently reliable. Here is a simple example of a divergent series and finding its "sum" with  evalf :

f:=unapply(`if`(n<10^10, 1/(2^n+3^n), 10^(-10)), n);

evalf(sum(f(n), n=1..infinity));

Obviously, Maple ignores the behavior of series for  n>=10^10 , and calculates the sum by using the first few terms.

For a rigorous solution of the problem can make the following two steps:

1) The proof of the convergence of the series easy to make, using  d'Alembert's ratio test .

2) For a rigorous calculation with a given accuracy is necessary to make an estimate of n-th remainder of the series R[n] = S - S[n] (here  S  is the sum of a series and  S[n]  is  n-th partial sum).

Step 1 we make with Maple:

u:=n->(n!/(2*n)!)^n:

limit(u(n+1)/u(n), n=infinity);

                                0

Step 2 by hand:

First we estimate  u[n+1]  and then estimate  R[n]  by an geometric series  sum(b*q^n, n=0..infinity)=b/(1-q)  with  q=1/2

Obviously,   

u(n+1)=((n+1)!/(2*(n+1))!)^(n+1)=(1/((n+2)*...*(2*n+2)))^(n+1)

<(1/(n+2))^(n^2+2*n+1)

This is the final estimate:   R[n] < 2/((n+2)^(n^2+2*n+1))

Example of use.

It is necessary to calculate the sum of the original series with an accuracy of 10^(-10) :

fsolve(2/((n+2)^(n^2+2*n+1))=10^(-10));

                        2.870597417

We see that in order to achieve the given accuracy is sufficient to find sum the first three terms of the series.

See the  result

In the assignment  

f:= x(t)^2 + y(t)^2:  

f  is an expression, not a function. This is inconvenient, if necessary to calculate the value of the function and its derivative at a point. See

f:= x(t)^2 + y(t)^2:

diff(f, t):

f(1);

                                   x(t)(1)^2+y(t)(1)^2

diff(f, t)(1);

                   2*x(t)(1)*diff(x(t),t)(1)+2*y(t)(1)*diff(y(t),t)(1)

It is better to define  f  and its derivative as  functions:

f:=t->x(t)^2+y(t)^2:

D(f);

                        t -> 2*x(t)*D(x)(t)+2*y(t)*D(y)(t)

f(1);

                                      x(1)^2+y(1)^2

D(f)(1);

                         2*x(1)*D(x)(1)+2*y(1)*D(y)(1)

Specific example:

f:=t^2;  # Assignment to the expression

                                      f := t^2 

f(1);  # Here Maple does not understand that it is necessary to substitute the number  1  instead of  t

                                       t(1)^2

diff(f, t);

                                         2*t

diff(f, t)(1);  # The same for the derivative

                                        2*t(1)

f:=t->t^2;  # The functional assignment

                                    f := t -> t^2

D(f);  # The derivative as a function

                                        t -> 2*t

f(1);  # Maple correctly calculates the value of the function  f  at the point  t=1

                                           1

D(f)(1);  # The same for the derivative

                                           2

 Edited.