All your constants need to have numerical values. So we need to know what those are.

The boundary conditions you give as:

bcs := y[1](0) = 0, y[2](0) = 1, y[3](0) = d[1], y[2](blt) = 0, y[3](blt) = 0, y[5](0) = 1, y[5](blt) = 0, y[4](0) = d[2], y[6](0) = d[3], y[4](blt) = d[4];

where blt is known, but d[1],..., d[4] are not. The six boundary conditions not containing the d's should determine the solution, so the values of y[3](0), y[4](0), y[6](0), y[4](blt) will be determined directly from the solution.
This has nothing to do with using a shooting method or not.

@emmantop Surely, you had a system of two odes of order 4 in f and order 2 in theta.
For some reason (which?) you want to turn this into a first order system using the variable names y[1],...,y[6].
Thus you must write down 6 first order equations of the form
diff(y[1](eta),eta)=y[2](eta);
diff(y[2](eta),eta)= whatever it is in terms of all the y[i]'s.
etc. up to and including
diff(y[6](eta),eta)= whatever it is in terms of all the y[i]'s

You have only two of these at the moment, but you have a line containing the translations from f,theta to the y[i]'s.
So the information from there has to be written in the form I showed above.

Basically, DEtools[convertsys] could have done that.
Still the question remains, why convert to a first order system? dsolve/numeric will do it itself (without telling you).

A suggestion if you insist on writing the original system of 2 equations (say eq1, eq2) in f and theta as a first order system:
Write down eq1 and eq2. Then do
DEtools[convertsys]({eq1,eq2},[],[f,theta],eta,Y,YP);
Then you get a list, whose first element is the first order system. Now write that in the form desribed above.

Before getting to the issue at hand you need to fix (at least) two things:

1. Missing multiplication sign in ode1 after
2. There is a theta(eta) appearing in ode2. That should probably be y[5](eta).

@surnizai The plot of the complex valued function x*y1 (defined in the complex x-plane,but restricted to the rectangle x=-2-I..2+I) is really a 3d-plot with the "x"-axis being Re(x), the "y"-axis being Im(x), and the "z"-axis being the absolute value of the function, i.e. abs(x*y1).

@Carl Love I tried the following:

res:=dsolve(eq);
ode:=op([2,2,1,1],res);
resb:=dsolve(eval(ode,epsilon=0));
eq3:=t=subs(_a=x,int(1/rhs(resb[1]),_a))+_C2;
eq4:=t=subs(_a=x,int(1/rhs(resb[2]),_a))+_C2;
sol1:=solve(eq3,x);
nops([sol1]);
simplify([sol1]); #Not simple
# It turns out that the other solution for b(a) leads to the same two solutions for x:
sol2:=solve(eq4,x);
nops([sol2]);
sol1[1]-sol2[1];
sol1[2]-sol2[2];

I don't know how to make it as simple as the direct answer from dsolve with epsilon=0.
But introducing two other arbitrary constants replacing _C1 and _C2 might get us somewhere closer.
indets(sol1[1],radical);
#This example is somewhat encouraging:
evalf(eval(sol1[1],{_C1=1,_C2=0}));

@Carl Love The chain rule is being used. x(t) is chosen as independent variable _a, _b(_a) is defined as diff(x(t),t).
Leaving out underscores then by the chain rule (and using epsilon=0 for simplicity):

db/da = (db/dt)*(dt/da) = x''(t)*(1/b(a)) = (x(t)-x(t)^3)*(1/b(a) .

Thus  b'(a)*b(a) - a + a^3 = 0. (*)

To get t we use dx/dt = b(a), thus using x(t)=a we get dt/da = 1/b(a) so that t = Int( 1/b(a),a) +C1.

The second constant C2 will appear when the ode (*) is solved.

Is it better than nothing? Good question, but at least the problem has been reduced to solving a first order ode and a subsequent integration. But t is then expessed in terms of a (i.e. x), t = t(x). So we have to invert to find x=x(t).
The stucture is described here:
??ODESOLStruc

@acer With the package LinearAlgebra both work:

restart;
eval(LinearAlgebra);
op(2,eval(LinearAlgebra));

I checked a few other packages, and found the same for

Statistics
Student
Units
LargeExpressions

There may be more.
####################### OK these are named modules ###################
To quote from the help page for named modules:
"Evaluating a named module definition causes, as a side effect, the assignment of the resulting module to the name, and the module name is stored in the attributes of each module member."

@Markiyan Hirnyk This seems to be endpoint problems caused by the fact that V(x,t) cannot be computed for values of x outside of the interval 0..1. Thus a two-sided limit cannot be computed by fdiff at x=0 or x=1.

V(1.00001,0.3);

Vx(0.99999,0.3);
fdiff(V,[1],[1,0.3]);
fdiff(V,[1],[0,.3]);
fdiff(V,[1],[0.0001,0.3]);
Vx(.0001,.3);

As it turns out your version has the same problem, since you exchanged the order of x and t in the definition of VVX.
I suppose you want
VVX := (b, a)-> fdiff(VV(x, a), x = b);
so that x is first as it is in V (and in your VV).

@northpole11 Here is a version, which is slightly different from Markiyan Hirnyk's.
Vx returns unevaluated for non-numeric input:

Vx:=proc(x,t) if not type([x,t],list(numeric)) then return 'procname(_passed)' end if;
   fdiff(V,[1],[x,t])
end proc;
Vx(.5,0.4567);
Vx(r,q);
plots:-animate(plot,[Vx(x,t),x=0..1],t=0..1);
plot3d(Vx(x,t),x=0..1,t=0..1);

@patient In the first few lines of your worksheet the definitions of eq1 and eqc has diff(w2,z) and diff(w2c,z) (respectively) on the right hand side. These will be evaluated to 0 and so eq1 and eqc will have zero right hand sides thereafter even if you later define w2 and w2c.
I guess that was not intended.
Solution: Replace w2 and w2c by w2(z) and w2c(z) in the equations defining eq1 and eqc.

I'm just guessing what you intend. The worksheet is without explanation and (at least to me) it is not obvious what you are trying to do.

@northpole11 In your worksheet example.mw the syntax for fdiff should be (example):

fdiff(V, [1], [.5, .4567]);
# or
fdiff(V(x, t), x, [x = .5, t = .4567]);

@tomleslie I'm sure I filed an SCR for the alpha-version present in Maple 18.02 probably in March.

@Rouben Rostamian  It is my experience after having tested quite a few examples from C.A.H. Paul, A Test Set of Functional Differential Equations (1994) that dsolve handles equations with delay dependent on y(t) quite well. That the plot always is wrong for t<t0 (start) I attribute to the implementation. Since the user knows that the history is constant (and what that constant is) this is not really a problem.
So maybe what is missing is a test for the delay becoming negative.

Although dsolve/numeric/delay assumes constant history (which we don't have here), we can easily get around that as will be seen here. The trick consists in differentiating the 'history' i.e. sin(t) and using a piecewise. Initially on an interval of length = Pi/2 (the delay) the equation is an ode, thus no history is required, just an initial condition. By the time we get to t=0 sufficient history has been built.

eq:=diff(y(t),t)=piecewise(t<0,cos(t),-y(t-Pi/2)):
res:=dsolve({eq,y(-Pi/2)=-1},numeric);
plots:-odeplot(res,[t,y(t)],-Pi/2..10);
plots:-odeplot(res,[t,y(t)-sin(t)],-Pi/2..10); #The maximal error is less than 1e-6.

I realize that you are using Maple 13, where solving delay differential equations numerically is not possible. It is in Maple 18 and Maple 2015.