In Maple versions Maple13 and later you could do:
op(op~([seq([filler[k],slf[k]],k=1..10)]));

In earlier versions (as well as later):
op(map(op,[seq([filler[k],slf[k]],k=1..10)]));

Example:
for i to 10 do slf[i]:=i; filler[i]:=" " end do:
printf("%s%d%s%d%s%d%s%d%s%d%s%d%s%d%s%d%s%d%s%d\n",op(map(op,[seq([filler[k],slf[k]],k=1..10)])));
# In Maple 2015.2 this prints:
1 2 3 4 5 6 7 8 9 10

I see that you have Maple 7. This version is not available to me at this time.

The formatting string could be done like this:
f:=cat("%s%d"$10);
printf(cat(f,"\n"),op(map(op,[seq([filler[k],slf[k]],k=1..10)])));

The elementwise operation ~ only works for so-called container objects, i.e. sets, lists, tables, arrays, vectors, and matrices.

A simple example made for the occasion just to show that map isn't at all obsolete:

u:=4^(1/2)+sqrt(x^2)/(1+x);
#u is not a container, but you may want to apply simplify to each individual term:
simplify(u) assuming x>0; #Applied to u
map(simplify,u) assuming x>0; # Applied termwise.

I believe that the elementwise operation was introduced in Maple 13.

Here is a numerical solution. I start by replacing x by y=x/L (and omega^2 by omega2):

restart;
Digits:=15; #Try later with Digits:=25 if you wish
L := 0.2e-5;
##
p1 := 868.13186813186813186*(diff(f1(x), x, x))+5.6428571428571428571*10^9*(diff(f2(x), x))+2.6043956043956043956*10^9*(diff(f3(x), x))-3.0384615384615384615*10^16*f1(x);
p2 := 303.84615384615384615*(diff(f2(x), x, x))-8.6813186813186813186*10^16*f2(x)-5.6428571428571428571*10^9*(diff(f1(x), x))-8.6813186813186813186*10^16*f3(x)+0.75000000000000000000e-5*omega^2*f2(x);
p3 := -7.2344322344322344322*10^(-17)*(diff(f3(x), x, x, x, x))+0.14468864468864468864e-1*(diff(f3(x), x, x))-2.6043956043956043956*10^9*(diff(f1(x), x))-8.6813186813186813186*10^16*f2(x)-8.6813910256410256409*10^16*f3(x)-6.2500000000000000000*10^(-25)*omega^2*(diff(f3(x), x, x))+0.75000625000000000000e-5*omega^2*f3(x);
##
bcs := {f1(0) = 0, f1(L) = 0, f2(0) = 0, f2(L) = 0, f3(0) = 0, f3(L) = 0, ((D@@1)(f3))(0) = 0, ((D@@1)(f3))(L) = 0}
##
sys:=subs(omega^2=omega2,{p1,p2,p3});
##
sys2:=PDEtools:-dchange({f1(x)=g1(y),f2(x)=g2(y),f3(x)=g3(y),x=L*y},sys,[g1,g2,g3,y]);
solve(sys2,{diff(g1(y),y,y),diff(g2(y),y,y),diff(g3(y),y$4)});
## Scaling omega2:
sys3:=subs(omega2=10^19*omega3,%);
## Boundary conditions:
bcs3:={g1(0) = 0, g1(1) = 0, g2(0) = 0, g2(1) = 0, g3(0) = 0, g3(1) = 0, D(g3)(0) = 0, D(g3)(1) = 0};
## Solution:
res3:=dsolve(sys3 union bcs3 union {(D@@2)(g3)(0)=.1},numeric,approxsoln=[omega3=1,g1(y)=y*(1-y),g2(y)=y*(1-y),g3(y)=y^2*(1-y)^2],maxmesh=1024,abserr=1e-5);
## Plots:
plots:-odeplot(res3,[seq([y,cat(g,i)(y)],i=1..3)],0..1);
## Finding the corresponding omega:
subs(res3(0),omega3);
#omega^2=omega2=10^19*omega3
omega=sqrt(10^19*%);

@iman Well, OK, you want an analytic solution.
You can use method=laplace on an initial value problem:
ics:={f1(0) = 0, f2(0) = 0, f3(0) = 0, D(f1)(0)=f1p,D(f2)(0)=f2p,D(f3)(0)=f3p,(D@D)(f3)(0)=f3pp,(D@@3)(f3)(0)=f3ppp};
where f1p, f2p, f3p, f3pp, and f3ppp are to be determined such that your bcs are satisfied.
But more than that: Since your problem has the trivial zero solution for any omega, you need to formulate a condition that ensures that nonzero solutions exist, which will only be for certain values of omega.

sol:=dsolve({p1,p2,p3} union ics,{f1(x),f2(x),f3(x)},method=laplace):

After that I say: Good Luck!

@Annonymouse Since nondimensionalization for a given physical system can be done in many ways, what kind of answer would you expect from such an automated approach?

Take as an example one I have also used recently on MaplePrimes.
This is a simple model for two competing species.

restart;
f:=(x,y)->x*(e1-s1*x-a1*y);
g:=(x,y)->y*(e2-s2*y-a2*x);
#The system:
sys:={diff(x(t),t)=f(x(t),y(t)),diff(y(t),t)=g(x(t),y(t))};
varchange:={x(t)=X*xi(tau),y(t)=Y*eta(tau),t=T*tau};
PDEtools:-dchange(varchange,sys,[tau,xi,eta]);
solve(%,{diff(xi(tau),tau),diff(eta(tau),tau)});
systemp:=expand(%);
#Now you want to pick X, Y, and T, so that the new system becomes dimensionless.
# This can be done in many ways.
# If you pick T=1/e1 then the coefficient of xi in the ode for xi will equal 1, but that isn't very good if you intend to include the possibility that e1=0 (or worse e1<=0).
# Proceeding, however, with the following choices:
#  T = 1/e1, X=e1/s1 and Y=e1/a1,
# which will assume that not only e1, but also s1 and a2 are nonzero, we get
eval(systemp,{T=1/e1,X=e1/s1,Y=e1/a1});
#We may give the coefficents in the ode for eta new names:
sys_nondim:=subs(a2=a*s1,s2=s*a1,e2=e*e1,%);
##############
## If you don't want to make these nonzero assumptions, then you could just say that T, X, and Y are positive constants with the same dimensions as t, x, and y, respectively.
# But then you won't be reducing the number of parameters in the system:
L:=subs(systemp,diff([xi,eta](tau),tau));
S:=map(coeffs,L,[xi(tau),eta(tau)])=~[seq(c[i],i=1..6)];
params:=indets(sys,name) minus {t};
S1:=solve(S,params);
newsys:=subs(S1,systemp);

restart;
u:=log[2*abs(x-a)](abs(x+a)+abs(x-a));
u1:=subs(abs(x-a)=b/2,abs(x+a)=c,u);
res1:=solve({u1 < 1,c>0}, c) assuming b>1;
res2:=solve({u1<1,c>0}, c) assuming b<1,b>0;
res1a:=eval(res1 union {b>1},{b=2*abs(x-a),c=abs(x+a)}); #Edited
sol1:=solve(res1a,{x,a});
res2a:=eval(res2 union {b<1,b>0},{b=2*abs(x-a),c=abs(x+a)}); #Edited
sol2:=solve(res2a,{x,a});
p1S:=seq(plots:-inequal(sol1[i],a=-1..1,x=-1..1),i=1..nops([sol1])):
p2S:=seq(plots:-inequal(sol2[i],a=-1..1,x=-1..1),i=1..nops([sol2])):
plots:-display(p1S,p2S);

You cannot use square brackets instead of parentheses.
Thus use:
eq2 := -7.500000000*10^7*(0.5e-2*(D(y))(t)-0.2500000000e-3*sin(t))/(0.5e-2*y(t)+0.2500000000e-3+0.2500000000e-3*cos(t))-1103.141678*(D(y))(t)*(1/(.5-y(t)))^1.4/(.5-y(t)) = 0.1e-1*((D@@3)(y))(t);

#Furthermore add maxfun=0 to your dsolve command:

sol2 := dsolve({eq2, ics}, numeric, maxfun = 0);

maxfun is the maximal number of function evaluations (evaluatens of the ode at different values of t).
Setting it at zero means actually setting it at infinity. You may want to set it instead at a high value like maxfun=10^8 so that the computation will terminate in finite time (and before we are all dead).
So try e.g.:
sol2 := dsolve({eq2, ics}, numeric, maxfun = 10^8);
plots[odeplot](sol2, 0 .. 10);

You have a missing semicolon after zahl+1.
Also I added to 10 in the beginning: adjust as you like, but I wanted this to end!
restart;
zahl:=1234567:
to 10 while sqrt(zahl)<> round(sqrt(zahl)) do  
   zahl:=zahl +1;  
     if sqrt(zahl)<> round(sqrt(zahl))then  
  print(zahl):  
 end if:
 end do;

You are using Euler's constant gamma as a variable name. That is as bad as using Pi as a variable name.
Try e.g. in your worksheet to do
evalf(bcs);
You will see that one of the boundary conditions has become .5772156649 = 1. , because evalf(gamma(0)) = evalf(gamma) = .5772156649.

The remedy is simple: Either from the very start use another name (g, e.g.) or when defining dsys do:
dsys := subs(gamma = g, {Eq1, Eq2, Eq3, Eq4, bcs});

First of all, try to remedy the immediate problem: that you must use a midpoint method. The error message says as much. Thus use either method=bvp[midrich] or method=bvp[middefer].
When you try again you get the error message
                Error, (in dsolve/numeric/BVPSolve) singularity encountered

If you isolate the highest derivatives in your system, i.e.
{diff(F(eta),eta$3),diff(E(eta),eta,eta),diff(K(eta),eta,eta),diff(theta(eta),eta,eta)}

you will see that K(eta)^2*E(eta) is the denominator in 3 equations and K(eta)^3*E(eta) in the fourth.

Since your boundary conditions include K(0) = 0, E(0) = 0, K(1) = 0, E(1) = 0 you may very well have an actual singularity at eta = 0 or eta=1 (or both).

Without assumptions it is difficult to see what can be done. Remember that Maple uses the principal square root, and that u a priori is complex.

Try these simplifications for u real:  u>1 or u<-1:

restart;
g := 1/2*u*(4*u-3)*sqrt(-(u-1)*(u+1)*(u^2-u-1))/sqrt(u*(u-1));
g1:=radnormal(g,rationalized);
combine(g1) assuming u>1;
combine(g1) assuming u<-1;
normal(%);
plot([Re,Im]~(g1),u=-2..2,-2..2,thickness=3); #Shows the necessity of splitting
solve(u^2-u-1,u);
evalf(%);

The two problems give different answers, but that is not surprising, since you have two quite different odes.

The first one is in fact:
fx:=int(exp(-(1/10)*t^2), t = 0 .. x);
dsolve({diff(y(x), x) = y(x)+fx, y(0) = 0}); # Symbolic solution for comparison with your numeric:
evalf(eval(%,x=1));
                    
The second is in fact:
f2:=eval(fx,x=1);
dsolve({diff(y(x), x) = y(x)+f2, y(0) = 0}); # Symbolic solution for comparison with your numeric:
evalf(eval(%,x=1));
                     
These results agree with your numeric results.

Only you can decide which of the two differential equations is relevant for you to solve.

Replace AA by:
AA2:=subs(Y=Y(r1,r2,r3),AA);

so the dependency on r1,r2,r3 becomes explicit. Then dchange works.

restart;
for i to 4 do
  v:=L/2*cos(w*t+phi[i]);
  Eq[i]:=u[i](t)=v
end do;

or even simpler

restart;
for i to 4 do   Eq[i]:=u[i](t)=L/2*cos(w*t+phi[i]) end do;

You can use select:

Pz:=select(s->vdot(s)>0,zeros);
##If you want the nonnegative ones too, you can do both at the same time:
Pz,Nz:=selectremove(s->vdot(s)>0,zeros);