restart;
sys:= diff(S(t), t) = -eta*K(t)*S(t)/(w*N*(S(t)+K(t))), diff(K(t), t) = eta*K(t)*S(t)/(w*N*(S(t)+K(t)))+S(t)*(-z*eta*alpha*K(t)^2+(-z*eta*alpha*S(t)-(eta*alpha^2*S(t)^2-2*N*C[max]*w*eta*alpha*K(t)+((-N*w+z)*alpha+N*C[max]^2*w*eta)*w*N)*upsilon)*K(t)+N*S(t)*w*alpha*upsilon*(N*w-z))/((K(t)^2*alpha*z+3*K(t)*S(t)*alpha*z+(2*S(t)*z*alpha+upsilon)*S(t))*w*N);
eq1:= -c2 + (K(t2)*S(t2)+w*N*S(t2))*z=0;
ics:=S(0)=100,K(0)=20;
params:=indets({sys,eq1},name) minus {t,t2};
##Now I set all parameters to 1 since I don't know what they are:
vals:=params=~1;
##Use your own of course.
eval([sys,eq1],vals);
eq2:=subs(t2=t,%[-1]);
sys2:=%%[1..2][];
#Because of my parameters eq2 never becomes satisfied so for illustration I used the equation
# lhs(eq2)=200 instead
res:=dsolve({sys2,ics,T(0)=0},numeric,discrete_variables=[T(t)::float],events=[[lhs(eq2)=200,T(t)=t]]);
res(5);
plots:-odeplot(res,[t,lhs(eq2)-200],0..0.1);

In Maple go to the menu Tools:  Tools/Tutors/Linear Algebra/Gaussian Elimination ...

Here is another version. It uses 'satisfies'.

createModule2 := proc(A::Matrix(square))
    local dim;
    dim := LinearAlgebra:-RowDimension(A);
    module()
        export det;
        det := (x::And(Matrix(square),satisfies(x->op([1,1],x)=dim))) -> LinearAlgebra:-Determinant(x)
    end module
end proc:

I shall just call your attention to two procedures CurveFitting:-PolynomialInterpolation and Statistics:-Fit.

data:=[[1030, 0], [380, 34], [270, 73], [240, 150], [85, 700], [22, 2000], [12, 5000]];
##This one uses collocation, i.e, the curve passes through all the points.
p:=CurveFitting:-PolynomialInterpolation(data, x);
plots:-pointplot(data,symbolsize=20); p1:=%:
plots:-display(plot(p,x=0..1030),p1); # Not great
## Best fit to exponential:
q:=Statistics:-Fit(a*exp(-b*x),data,x);
plots:-display(plot(q,x=0..1030),p1); #Not particularly great either

So think of some reasonable model and use Statistics:-Fit.

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(%);