The second time around z is no more a symbol.
I understand this to be a small version of something more complicated. Is there a need to introduce symbols in Z?
###
If not then of course there is no problem:
restart;
for k from 1 to 2 do Z:=Matrix(2,2); Z[1,1]:=1; print(Z); end do;
###
Anyway to solve the problem as presented you could do either one of these:

restart;
for k from 1 to 2 do Z:=Matrix(2,2,symbol=z); z[1,1]:=1; print(Z); z:='z'; end do;
###
restart;
for k from 1 to 2 do Z:=Matrix(2,2,symbol=z); Z[1,1]:=1; print(Z); end do;

Do you actually have a nonlinear eigenvalue problem, where unfortunately for dsolve you have already found an eigenvalue and that is inserted into your system after which there is no parameter only concrete values for the constants.
If so it is better to let Maple find an eigenvalue and a corresponding solution at the same time.

Let me start by giving a very simple example after which I shall make an experiment with your problem.
restart;
ode:=diff(x(t),t,t)+lambda*x(t)^2=0; #Nonlinear ode
res:=dsolve({ode,x(0)=0,x(1)=0,D(x)(0)=1},numeric); #Extra requirement D(x)(0)=1 <> 0
res(0); #You see the eigenvalue
plots:-odeplot(res,[t,x(t)],0..1); #Corresponding solution for x
#Now insert the eigenvalue found into ode:
ode1:=subs(res(0)[-1],ode);
#Try solving with homogeneous boundary conditions:
res1:=dsolve({ode1,x(0)=0,x(1)=0},numeric);
res1(0.5); #Indicates clearly enough already that you just have the zero solution
res2:=dsolve({ode1,x(0)=0,D(x)(0)=1},numeric); #Using the extra condition instead of x(1)=0.
plots:-odeplot(res2,[t,x(t)],0..1);
################################
Now let us try an experiment with your concrete system containing no parameter lambda.
We shall rather arbitrarily insert a parameter lambda into your system and shall add an extra (nonzero) boundary condition.

I use the same ideas I used in my other answer, but I repeat it here.

restart;
dsys3 := {0.1873004296e-2*(diff(f3(x), x, x, x, x))-.3714109950*(diff(f3(x), x, x))+13.96314014*(diff(f1(x), x))+1.822500000*10^(-10)*(diff(f2(x), x, x))-4.888800000*10^(-7)*f2(x)+22.72337812*f3(x) = 0, -901.7765285*(diff(f1(x), x, x))+3.153274902*10^(-10)*(diff(f2(x), x, x, x))+0.4807756008e-5*(diff(f2(x), x))-27.92641797*(diff(f3(x), x))+18074.25057*f1(x)-3.100642108*(diff(f3(x), x, x))*(diff(f3(x), x))+7.74301535*(diff(f3(x), x))*f3(x) = 0, 0.2249991345e-1*(diff(f2(x), x, x, x, x))-227.1861332*(diff(f2(x), x, x))+29537.78948*f2(x)-2.841067252*10^(-11)*(diff(f1(x), x, x, x))-3.392672540*10^(-7)*(diff(f1(x), x))+7.650000000*10^(-11)*(diff(f3(x), x, x))-1.400000000*10^(-7)*f3(x)-5.500000000*10^(-7)*f3(x)^2+9.278333333*10^(-13)*(diff(f3(x), x))^2+3.300000000*10^(-13)*(diff(f3(x), x, x))*f3(x) = 0, f1(0) = 0, f1(1) = 0, f2(0) = 0, f2(1) = 0, f3(0) = 0, f3(1) = 0, (D(f2))(0) = 0, (D(f2))(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0};
for i to 3 do dsys3[i] end do;
diff(dsys3[2],x);
eq:=subs(solve(dsys3[3],{diff(f2(x),x$4)}),%);
sys:=solve({dsys3[1],dsys3[3],eq},{diff(f1(x),x$3),diff(f2(x),x$4),diff(f3(x),x$4)});
dsys3[4..-1];
nops(%);
cdn:=eval(convert(dsys3[2],D),x=0);
#The insertion of lambda will be done here (obviously this is completely arbitrary and just an illustration):
sysE:={sys[1],sys[2],lhs(sys[3])=lambda*rhs(sys[3])};
resE:=dsolve(sysE union dsys3[4..-1] union {cdn, D(f1)(0)=1},numeric);
resE(0);
#Eigenvalue is -.493863843115919 i.e. not 1 so we don't have your original system.
#However, I inserted this lambda arbitrarily without knowing your original problem, so lambda=1 was hardly to be expected.
plots:-odeplot(resE,[[x,f1(x)],[x,f2(x)],[x,f3(x)]],0..1,thickness=3);

Now that we have the constants, we can begin by observing that solving for the highest derivatives isn't possible.

solve(dsys3[1..3],{diff(f1(x),x$3),diff(f2(x),x$4),diff(f3(x),x$4)});
##However by manipulating the system some, e.g. as follows you can solve the resulting system for the highest derivatives.
Please pay attention to if the ordering in the set dsys3[1..3] is such that dsys3[2] is the ode without any 4th derivative. It is assumed in the following:
Since we are differentiating dsys3[2] there will be the constraint dsys3[2] to keep in mind.
restart;
dsys3 := {0.1873004296e-2*(diff(f3(x), x, x, x, x))-.3714109950*(diff(f3(x), x, x))+13.96314014*(diff(f1(x), x))+1.822500000*10^(-10)*(diff(f2(x), x, x))-4.888800000*10^(-7)*f2(x)+22.72337812*f3(x) = 0, -901.7765285*(diff(f1(x), x, x))+3.153274902*10^(-10)*(diff(f2(x), x, x, x))+0.4807756008e-5*(diff(f2(x), x))-27.92641797*(diff(f3(x), x))+18074.25057*f1(x)-3.100642108*(diff(f3(x), x, x))*(diff(f3(x), x))+7.74301535*(diff(f3(x), x))*f3(x) = 0, 0.2249991345e-1*(diff(f2(x), x, x, x, x))-227.1861332*(diff(f2(x), x, x))+29537.78948*f2(x)-2.841067252*10^(-11)*(diff(f1(x), x, x, x))-3.392672540*10^(-7)*(diff(f1(x), x))+7.650000000*10^(-11)*(diff(f3(x), x, x))-1.400000000*10^(-7)*f3(x)-5.500000000*10^(-7)*f3(x)^2+9.278333333*10^(-13)*(diff(f3(x), x))^2+3.300000000*10^(-13)*(diff(f3(x), x, x))*f3(x) = 0, f1(0) = 0, f1(1) = 0, f2(0) = 0, f2(1) = 0, f3(0) = 0, f3(1) = 0, (D(f2))(0) = 0, (D(f2))(1) = 0, (D(f3))(0) = 0, (D(f3))(1) = 0};
for i to 3 do dsys3[i] end do; #JUST FOR HAVING A LOOK at the equations
diff(dsys3[2],x);
eq:=subs(solve(dsys3[3],{diff(f2(x),x$4)}),%);
sys:=solve({dsys3[1],dsys3[3],eq},{diff(f1(x),x$3),diff(f2(x),x$4),diff(f3(x),x$4)});
dsys3[4..-1]; #The boundary conditions
nops(%);
cdn:=eval(convert(dsys3[2],D),x=0); #The above mentioned constraint
res:=dsolve(sys union dsys3[4..-1] union {cdn},numeric);
res(0.5);
plots:-odeplot(res,[x,f3(x)],0..1);
#Now belatedly we notice that there is a trivial solution:
odetest({f1(x)=0,f2(x)=0,f3(x)=0},dsys3[1..3]);
############
Thus you are I assume looking for the possibility of a nontrivial solution?

I think the small parameter method is something like the following, where because of only two new independent variables x1 and x2 we go to second order in epsilon. It is important of course to notice that x(t) = 1 is a solution (in particular when epsilon=0).

restart;
eq:= diff(x(t),t$2) - epsilon*(alpha*x(t)^4 - beta)*diff(x(t),t) - x(t) + x(t)^3;
eq2:= eval(eq, x(t)= 1 + epsilon*x[1](t) + epsilon^2*x[2](t));
collect(eq2,epsilon,factor);
ode1:=coeff(%,epsilon);
ode2:=coeff(%%,epsilon^2);
res:=dsolve({ode1,ode2},explicit);
#To shorten:
combine(res);

#################Numerical comparison for epsilon = 0.1 ##########
res1:=dsolve({eval(eq,{alpha=1,beta=1,epsilon=.1}),x(0)=1.11,D(x)(0)=0},numeric); #1.11 = 1+eps+eps^2
plots:-odeplot(res1,[t,x(t)],0..15); p1:=%:
##Using numerical solution here for convenience in the comparison
res2:=dsolve({ode1,eval(ode2,{alpha=1,beta=1}),x[1](0)=1,x[2](0)=1,D(x[1])(0)=0,D(x[2])(0)=0},numeric);
plots:-odeplot(res2,[t,1+0.1*x[1](t)+0.01*x[2](t)],0..15,color=blue); p2:=%:
plots:-display(p1,p2);
##To use the symbolic solution for x1,x2 just remove 'numeric':
res2a:=dsolve({ode1,eval(ode2,{alpha=1,beta=1}),x[1](0)=1,x[2](0)=1,D(x[1])(0)=0,D(x[2])(0)=0});
X:=subs(res2a,1+0.1*x[1](t)+0.01*x[2](t));
p3:=plot(X,t=0..15,color=blue);
plots:-display(p1,p3);

Was your intention with p something like this:

simplify(piecewise(abs(z)>6^(2/3)/4,piecewise(z<0, eq30, eq31),piecewise(z<0, eqc20, eqc21)));
p:=unapply(%,z,t);

eq:=diff(y(x),x)=cos(x)+y(y(x)-2);
res:=dsolve({eq,y(0)=1},numeric,delaymax=1);
plots:-odeplot(res,[x,y(x)],0..1);

#If you try res(2) you get an error message saying:
res(2);
Error, (in res) cannot evaluate the solution further right of 1.0058110, probably a singularity

You could use Spline from the CurveFitting package:

res:=CurveFitting:-Spline(Data1,x);
m,M:=(min,max)(map2(op,1,Data1));
plot(Data1,style=point,symbolsize=20,color=blue); p1:=%:
plot(res,x=m..M); p2:=%:
plots:-display(p1,p2);

Somehow `dsolve/numeric/bvp/convertsys` fails to isolate YP[4] when converting to a first order system. I shall show that below. YP[4] stands for diff(f(eta),eta$4).
A workaround is to do this isolation before we try dsolve.

Thus e.g. one could do:

ode12:=solve({ode1,ode2},{diff(f(eta),eta$4),diff(theta(eta),eta,eta)});
sys:=eval(ode12 union {bcs},E);
res:=dsolve(eval(sys,epsilon=L[1]),numeric,method=bvp[midrich],maxmesh=1024);
plots:-odeplot(res,[eta,f(eta)],0..5);

##### Another workaround (simpler):
sys:=convert(sys,rational);
## A third workaround:
sys:=collect(sys,diff);

####################
To exhibit the problem you can begin with the redefined 'sys' from above:
`dsolve/numeric/bvp/convertsys`(eval(sys[1..2],epsilon=L[1]),sys[3..-1],[f(eta),theta(eta)],eta,cont);
#Obviously here YP[4] is isolated on the left, since it started out that way.
# Now try the nonisolated system:
sys:=eval([ode1, ode2, bcs], E);
`dsolve/numeric/bvp/convertsys`(eval(sys[1..2],epsilon=L[1]),sys[3..-1],[f(eta),theta(eta)],eta,cont);
## You will notice that for i=1,2,3,5,and 6, YP[i] has been isolated, but an equation defining YP[4] implicitly has been left unsolved. I don't believe that that is intentional: I tried changing f(0)=0 to f(0)=1 and that didn't make any difference to convertsys.
So my guess is that it is indeed a bug in `dsolve/numeric/bvp/convertsys`.
####
It appears that the problem is in line 55 of `dsolve/numeric/bvp/convertsys` where a (global) procedure `dsn/solve` fails to isolate YP[4].
Finally and ironically, the problem in `dsn/solve` is with the isolate command in line 30 where isolate cannot isolate YP[4].

Under the assumption that b^2-4*c < 0 clearly the integral must be real (except of course for an additive constant: the integral is indefinite).
As an example I plot the value for b=2, c=2, where b^2-4*c = -4 < 0:
res:=int((y^4+b*y^2+c)^(-1/2),y);
plot(eval(res,{b=2,c=2}),y=-8..8);

Notice that the second derivatives of alpha and beta only appear in eq2. Thus solving (algebraically) for the second derivatives cannot be done as in

solve({eq1,eq2}, diff({alpha(t),beta(t)},t,t));

By using the word 'algebraically' I mean that if we consider the derivatives of alpha and beta of orders 0,1,2 as just six variables (like a0,a1,a2,b0,b1,b2) we cannot solve for {a2,b2} in terms of {a0,a1,b0,b1}.
We can get around it though: See the second method below.

Thus in order to find expressions for the second derivatives we could differentiate eq1 or let Maple do it. However, in differentiating eq1 information is lost as any constant on the right becomes a zero. Thus this equation must be kept at least evaluated at one point.

Anyway, Maple can do this itself like this:

solution := dsolve([eq1, eq2, CI], numeric, range = 0 .. 2,method=rkf45_dae);
#You will get an error about the initial conditions not satisfying the constraint, and in fact they don't:

eval(convert(eq1,D),t=0);
eval(%,{CI});

Thus you must clear up that problem, i.e. change the initial conditions.
################################ Alternative solution ################
Alternatively you can do like the following. Here we avoid any differentiation.
Notice that no requirements can be put on the derivative of beta at zero. This reflects the constraint problem above.
eq_beta:=solve(eq1,{diff(beta(t),t)});
eq3:=eval(eq2,%);
solution:=dsolve({eq3,op(eq_beta),alpha(0) = 0, beta(0) = 0,D(alpha)(0) = w[0]},numeric);
plots:-odeplot(solution,[[t,alpha(t)],[t,beta(t)]],0..0.9);

The _Z in the RootOf result is just a dummy variable, the unknown in the expression inside RootOf. The meaning is that the result is a root in that expression, but Maple either could not express it analytically or didn't 'want' to unless forced to as in the case of the 4 solutions to a quartic polynomial (see allvalues).
In your case there is no chance for an analytical solution (i.e. a  formula for the solution).

If you set

res:=solve(equ,c);

then you can use evalf to get a numerical evaluation of (one of) the root(s):
evalf(res);

You will find a result that looks like it is likely that c=3100 is an exact solution. It is clear from your equation that c=3100 is indeed a solution, since 3100^2 = 9610000, so c=31000 makes both sides of the equation zero.

This solution is apparently not the one you want, so you could try fsolve with a start guess for c:

fsolve(equ,c=13000);
You will find a c-value close to zero. Not the one you want either.

For the square roots all to be real you need  c > 6300 thus it makes sense to simplify equ to:

equa:=simplify(equ) assuming c>6300;

You could try plotting the difference between the two sides (of equ or equa):
plot((rhs-lhs)(equa),c=6300..15000);

After that you should ask yourself if there are indeed any solutions or if you just wrote down the equation incorrectly.

The operator Mf of multiplication by the function f can be defined like this:

Mf:=g->f*g;
#Applying Mf to a function h (in the mathematical sense, not Maple)
(Mf(h))(x); #The parenthesis around Mf(h) i.e. (Mf(h)) are superfluous, but are added for clarity of meaning.
((D-Mf)@@2)(g)(x);
expand(%);

#An operator M which maps functions onto multiplication operators:

M:=f->(g->f*g);
#Thus M(f) is the operator which acts like this:
M(h)(g)(x);
((D-M(f))@@2)(g)(x);
expand(%);

########Multiplying by a constant function:
M(89)(g)(x);
a:=x->k; #The function a defined by a(x)=k for all x
M(a)(g)(x);

############## Followup question
#To accomplish your intended result you can use M above:

restart;
M:=f->(g->f*g);
id:=(x,m)->x;
M(id)(g)(y,n); #test
(M(id)@D[1])(g)(x,m); #test
((M(id)@D[1])@@2)(F)(x,m);

#Alternatively use id but not M:
Mx:=g->id*g;
((Mx@D[1])@@2)(F)(x,m);

You can use approxsoln (after correcting the B-problem pointed out by Carl):

AP:=NULL:
for k to 6 do
R := dsolve(eval({Eq1, Eq2, bcs1, bcs2}, B = L[k]), [f(eta), theta(eta)], numeric, output = listprocedure,maxmesh=512,AP);
AP:=approxsoln=R;
X1 || k := rhs(R[3]); X2 || k := rhs(R[4]); Y1 || k := rhs(R[5]); Y2 || k := -rhs(R[6])
end do;

The print command at the end then gives me the values:

[-1.41421356296431, -1.43879570677334, -1.46322186668648, -1.48747814527757, -1.51155337283017, -1.55912751987572]

Incidentally, there is no reason for the 'print'. Just writing
[(X2 || (1 .. 6))(0)];

will do.

I'm answering your question about tangency points.

Since fsolve is rather reluctant to return a result if it not perfect, or if it doesn't appear so to fsolve, I once wrote a procedure to pick the best solution found by fsolve in a failed attempt.
It is given below and called BestPoint.
It requires as input either a procedure or a list of procedures with option remember.
It finds the input argument(s) corresponding to the least value of an error function applied to the output of the procedure(s).
To test it just copy and paste the whole thing and execute it in the order given.
I'm applying it to your tangency problem.
The procedure on which BestPoint will be used is called q below.

restart;
BestPoint:=proc(q::{procedure,list(procedure)},{errorfunction::procedure:='default',procedurelist::truefalse:=false,
number::posint:=`if`(not procedurelist and q::list(procedure),nops(q),1)}) local Q,f,x,i,T,T1,E,R,S,inf;
 inf:=proc(eq) lhs(eq)=evalindets~(rhs(eq),Not(numeric),x->infinity) end proc;
 if errorfunction ='default' then f:=unapply(add(x[i]^2,i=1..number),seq(x[i],i=1..number)) else f:=errorfunction end if;
 if not procedurelist then
    if not type(q,list) then Q:=[q] else Q:=q end if;
    T:=op~(4,eval~(Q)); #Remember tables
    if T=[] then error "No remember tables" end if;
    T:=table~(map(inf~,op~(eval~(T)))); #replacing nonnumeric entries with infinity
    E:=[entries]~(T,nolist);
    R:=ListTools:-Enumerate(f~(op(E)));
    T1:=T[1]    
 else
    T:=op(4,eval(q));
    if T=NULL then error "No remember table" end if;
    T:=table(map(inf,op(eval(T)))); #replacing nonnumeric entries with infinity
    E:=[entries](T,nolist);
    if number<>nops(E[1]) then error "Wrong or missing argument: number = n" end if;
    R:=ListTools:-Enumerate((f@op)~(E));
    T1:=T
 end if;
 S:=sort(R,(x,y)->abs(x[2])<abs(y[2]));
 print(ErrorFunction=eval(f),Error=S[1,2]);
 indices(T1)[op([1,1],S)]
end proc:

#Your tangency problem.
Output from q is a list of two real numbers. These numbers we want to make as small as possible and will use fsolve on q in order to do that. This will fail in many cases, but then we use BestPoint.
On my run of the H-loop at the bottom (with Digits=10) fsolve succeeded 5 times out of 14 in Maple 2015.

q := proc (x,y,H) option remember;
local f0,f1, f2, k;
  f0:= 0.585=abs(k*(0.585+H*cos(20*Pi/180))+0.585-H*sin(20*Pi/180))/sqrt(k^2+1);
  f1:=fsolve(f0,k)*x-y;
  f2:= (x-0.585-H*cos(20*Pi/180))^2+(y+0.585-H*sin(20*Pi/180))^2-0.585^2;
  evalf([f1,f2])
end proc;
q1:=proc(x,y) q(x,y,H)[1] end proc;
q2:=proc(x,y) q(x,y,H)[2] end proc;
#fsolve cannot take q as input; that is the basic reason for q1,q2.

H:=-0.18: sol := fsolve([q1,q2]);
res:=table();
res[1]:=BestPoint(q,procedurelist,number=2);
printlevel:=2:
i:=1:
for H from -0.18 by 0.03 to 0.18 do
  i:=i+1;
  sol:=fsolve([q1,q2],res[i-1][1..2]);
  if not type(sol,list(numeric)) then res[i]:=BestPoint(q,procedurelist,number=2) else res[i]:=[op(sol),H] end if;
  forget(q)
end do;
#The results:
seq(res[i],i=1..14);
#Checking results:
seq(q(op(res[i])),i=1..14);

There is no such thing as the particular solution. Any solution is a particular solution. Thus the simplest way to get a particular solution without arbitrary constants is to use dsolve after which one can set the arbitrary constants to whatever one wishes. I shall choose to set both to 0.

restart;
f := diff(u[1](t), t, t)+u[1](t)+(3/4)*a^3*cos(beta+t)+(1/4)*a^3*cos(3*beta+3*t)=0;
dsolve(f);
eval(%,{_C1=0,_C2=0});