I do not see any bug. Just Maple not fully solve the problem. Here is a workaround:

restart;

Sol:=rsolve({f(n)=0.5*f(n-1)+0.5*f(n+1), f(0)=1, f(6)=0},f(n));

solve(eval(%,n=0)=1, f(5));

f:=unapply(eval(Sol, f(5)=%), n);

Increase the accuracy of the calculations. For example, if  Digits:=50  everything is all right .

Addition: it's interesting that if you use the  solve  command instead of  Eigenvectors  command, you get the best results:

Digits := 50:
p := unapply(CharacteristicPolynomial(M, x), x):
E := [solve(p(x))];
for k to 8 do p(E[k]) end do;

Obviously, this set is the boundary of the intersection  three regions in the space  -x+y+z <= 1, x-y+z <= 1, and x+y-z <= 1,  that is, the boundary of a polyhedral angle.

A similar example in the plane - the specification of an equilateral triangle by one equation:

plots[implicitplot](max(-y, y-sqrt(3)+sqrt(3)*x, y-sqrt(3)-sqrt(3)*x)=0, x=-1.5..1.5, y=-1..2, thickness=3, gridrefine=3, scaling=constrained);

As for relatively close issues,  see this thread. 

Here is an other solution

ContoursWithLabels(exp(2*x/(x^2+y^2)), -2 .. 2, -2 .. 2, {1, 2, 10, 0.1, 0.5}, [axes = box, color = red, thickness = 2]);

 Edited.  The original function has a singularity at the origin. If you "cut out" it, you can get a pretty good coloring:

f := `if`(x^2+y^2 < 0.015, undefined, exp(2*x/(x^2+y^2))):

ContoursWithLabels(f, -2 .. 2, -2 .. 2, {1, 2, 10, .1, .5}, [axes = box, color = "DarkBlue", thickness = 2], Coloring = [colorstyle = HUE, colorscheme = ["White", "Red"], style = surface]);

restart;

assign(seq(Equ[i]=(u[i](t)=L/2*cos(w*t+phi[i])), i=1..4));

For example:

Equ[1];

                      u[1](t) = 1/2*L*cos(t*w+phi[1])

and so on.

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)