I tried this and it worked:

restart;
fd:=["F:/file1.m", "F:/file2.m"];

for i from 1 to nops(fd)
do
cat(var,i):=12345*i;
save cat(var,i),fd[i];
end do;

restart;
fd:=["F:/file1.m", "F:/file2.m"];
for i from 1 to nops(fd)
do
read fd[i];
end do;
seq(cat(var,i),i=1..nops(fd)): %;

Change BCs := [x(0), x(0.2e-1)-0.15e-2] to

BCs := x(0), x(0.2e-1)-0.15e-2;

But you still have problems, but they are not syntactical anymore.

method=bvp[midrich] could be tried, but now you have other problems.

Try

sol := dsolve({BCs, eq}, {x(t)}, type = numeric,method=bvp[midrich],initmesh=1024);

plots:-odeplot(sol,[t,x(t)],0..0.02);

In the standard interface you could use Vector ( = Vector[column]):

Vector([x^2-1,x^3+1]):
factor~(%);

Here is an attempt that looks like it is working. Maybe there is a better way.

I begin by changing variables. The change x= z*d/2 is just for convenience.

The idea behind the other change y(x)=u(z)*y0 is that y0 should be considered y(0) later. Your problem as stated obviously has the trivial solution y(x) = 0 for all x. You want to find nontrivial solutions. Therefore that change.

The problem then is: Given d find y0 so u(0)=1.

restart;
ode:= diff(y(x),x$2) + 1/x *diff(y(x),x) + 0.6*y(x) + y(x)^3 - y(x)^5  = 0;
PDEtools:-dchange({x=z*d/2,y(x)=u(z)*y0},ode,[z,u]);
ode2:=expand(%*d^2/4/y0);
#New boundary conditions
bcs2:=D(u)(0) = 0,   u(1) = 0.8;
#Don't forget option remember in the following procedure, which outputs u, so if you want y (as a function of z) then you must multiply by y0.
pu:=proc(dd,y00) option remember; local ode3,U,L;
 if not type([dd,y00],list(numeric)) then return 'procname(dd,y00)' end if;
 ode3:=eval(ode2,{d=dd,y0=y00});
 L:=dsolve({ode3,bcs2},numeric,method=bvp[midrich],output=listprocedure);
 subs(L,u(z))
 end proc:

plot(pu(80,1.19),0..1);
seq(pu(80,1.192+c/10000)(0),c=0..9);
plot(1.1925*pu(80,1.1925),0..1);
#I didn't succeed in getting fsolve to produce anything, but see implicitplot below.
fsolve(pu(80,y0)(0)=1,y0=1.9);
plots:-animate(plot,[pu(80,y0),0..1,thickness=2],y0=1.1..1.3,frames=10);
plots:-implicitplot(pu(d,y0)(0)=1,d=2..100,y0=.8..1.2);

#Added after lunch (extracting y0-d-information from implicitplot and using that y0 is constant eventually):

imp:=plots:-implicitplot(pu(d,y0)(0)=1,d=2..10,y0=.8..1.2): imp;
A:=op(indets(imp,Array)):
A[99,1..2];
Y0:=unapply(piecewise(d<=10,CurveFitting:-Spline(A,d),A[99,2]),d):
plot(Y0(d),d=2..100);
plot([seq(Y0(k)*pu(k,Y0(k)),k=[2.1,2.2,2.5,3,5,10,100])],0..1,view=0..1.2);

It seems to turn up when using 2d-input. When using 1d-input ("Maple Notation") it doesn't show up.

Since the other replies use numerical solution of the ode, I thought I would point out that of course you can plot the exact solutions, when they are available, which admittedly often is not the case. But certainly in the example you mention. To make the plots more exciting I have used animation.

with(plots):
de := diff(y(x),x,x) + 4*y(x) = sin(2*x);
res:=dsolve({de,y(0)=y0,D(y)(0)=y1});
animate(plot,[eval(rhs(res),y0=0),x=0..3*Pi,caption=typeset([y(0)=0,D(y)(0)=y1])],y1=-2..2,trace=24,paraminfo=false);

#And now an animation inside another animation:
animate(animate,[plot,[rhs(res),x=0..3*Pi,caption=typeset([y(0)=y0,D(y)(0)=y1])],y0=-2..2,frames=10,paraminfo=false],y1=-2..2,frames=10,paraminfo=false);

I'm pretty sure that the problem is your use of !!! (Execute the whole worksheet).

Remember that it is the last assignment (timewise) that counts. If after executing the whole worksheet you return to a line above, then Maple uses that last assignment which may not correspond to the situation in that line.

Executing the whole worksheet is very often not a good idea, unless you have no intention to return to a previous line. 

If you do

dsolve(ode);

Maple will find the same solutions you could find by using the usual trick of multiplying the ode by y' and integrating to obtain

(*) diff(y(x),x) = (+/-) sqrt(a*y(x)^2-1/2*y(x)^4+1/3*y(x)^6+C)

Using the boundary condition D(y)(0) = 0 gives you a connection between C and y(0). Using y(2) = 0 when separating variables in (*) gives you an equation between y(0) and a by integrating from x = 0 to x = 2, something like

int(1/sqrt(a*y^2-y^4/2+y^6/3-(a*y0^2-y0^4/2+y0^6/3)),y=y0..0)=2; (using the plus sign from (*)).

Now you have to do integrate (!!!) and solve for y0 in terms of a. That is most likely not possible.

Even solving for y0 numerically when given a seems unlikely to succeed.

Do you have any reason to believe that there are nontrivial solutions to your boundary value problem for some value of a (or all)?

Looking into the numerical solutions it seems that there are solutions for negative values of a, e.g. taking a = -0.5 then there are solutions for y0 near -0.43 and -1.

Doing

dsolve({eval(ode,a=-0.5),D(y)(0) = 0, y(0) = -1},numeric);
plots:-odeplot(%,[x,y(x)],0..2);

seems to confirm that.

To answer your last question:

Solving D(L)(z) =( L(z)/(1+z) )+( (1+z)/H ), with H=Y(z) being a known numerically determined function, you can use the option known = Y.

You first have to make sure that Y returns unevaluated when not receiving numeric input (that is  'procname(z)' below).

restart;
with(plots):
eq := z-> (H^2+(-1)*.27*(1+z)^3-(1/20000)*(1+z)^4)/(H^2)^.1 = .7299500000;
Y := z->if not type(z,numeric) then 'procname(z)' else fsolve(eq(z), H=1) end if:
p := dsolve({D(L)(z) = L(z)/(1+z)+(1+z)/Y(z),L(0)=0}, type=numeric, range=0..5,known=Y);
odeplot(p,[z,L(z)],0..5);

It might be better or at least as good to take the previously mentioned differential equations approach and now solving a system of two ode's numerically:

yp := implicitdiff(eq(z), H, z);
ode:=diff(H(z),z)=subs(H=H(z),yp);
p := dsolve({D(L)(z) = L(z)/(1+z)+(1+z)/H(z), ode,L(0)=0,H(0)=1}, type=numeric, range=0..5);
odeplot(p,[z,L(z)],0..5);

I looked at your worksheet, but didn't run it.

However, just before your problems start you define C to be IdentityMatrix(2). This is very likely to cause all kinds of problems because C has already been used as C[alpha.f].

It is never a good idea to use a variable name and an indexed version of it independently.

Try this:

restart;
allvars:={C[alpha.f]=20,C[alpha.r]=24}:
C:=LinearAlgebra:-IdentityMatrix(2):
allvars;
Error, bad index into Matrix

Why start another thread and make it difficult for people to follow?

Anyway, the syntax is

implicitplot(eq(x), x=0..10,y=0..120);

and gives a curve above the x-axis. I am assuming that, as previously,

eq:=x->y^(-1/5)*(y^2-(1+x)^3-(1+x)^4)=1;

Using the example from one of your other questions you can use either fdiff or implicitdiff to plot anything with diff(y(x),x) in it.

To plot y(x) + 1 is straightforward.

eq:=x->y^(-1/5)*(y^2-(1+x)^3-(1+x)^4)=1;
Y:=x->fsolve(eq(x),y);
plot(Y,-10..10);
plot(Y+1,-10..10);
plot(fdiff(Y,[1],x),x=-2..2);
yp:=implicitdiff(eq(x),y,x);
plot(eval(yp,y='Y(x)'),x=-2..2);

Rewriting for convenience:

eq:=x->y^(-1/5)*(y^2-(1+x)^3-(1+x)^4)=1;

The numerical solution procedure using fsolve:

Y:=x->fsolve(eq(x),y);
plot(Y,-10..10);

Did you try fsolve?

fsolve({fff},{T, x[22], x[23], x[31], x[32], y[21]});

With a starting point (maybe you expect T > 0?)

fsolve({fff},{T=12, x[22], x[23], x[31], x[32], y[21]});

u1:=sin(sqrt(K)*t/sqrt(N));
subs(K=N*omega^2,u1);
simplify(%) assuming positive;

#No problem with the next version

u2:=sin(sqrt(K/N)*t);
subs(sqrt(K/N)=omega,u2);

The necessity of assumptions is due to the fact that sqrt(K/N) is not always equal to sqrt(K)/sqrt(N).

As an example consider K=1 and N = -1, in which case sqrt(K/N) evaluates to i, whereas sqrt(K)/sqrt(N) evaluates to 1/i = -i.