@9009134 I believe that you can trust the numerical result. The method itself is perfectly legitimate (recall the simple (linear) example I gave above).
I would do a simple test for consistency, by which I mean: If I feed RES as an approximate solution to dsolve/numeric/bvp is the solution recovered?
It is:

plots:-odeplot(RES,[[x,1000*f1(x)],[x,1000*f2(x)],[x,f3(x)]],0..1); p1:=%:
#omega = 948.876522424828
sol:=dsolve(dsys3 union {(D@@4)(f3)(0)=-1}, numeric,approxsoln=RES);
plots:-odeplot(sol,[[x,1000*f1(x)],[x,1000*f2(x)],[x,f3(x)]],0..1); p2:=%:
sol(0);
#omega = 948.876522440829
plots:-display(p1,p2); #Indistinguishable

@9009134 Just copy and paste the following code into a fresh worksheet:

restart;
dsys3 := {15*(diff(f3(x), x, x, x, x, x, x))+16*(diff(f3(x), x, x, x, x))+17*(diff(f3(x), x, x, x, x))+18*(diff(f3(x), x, x))+19*(diff(f3(x), x, x))+20*(diff(f3(x), x, x))+21*(diff(f1(x), x, x, x))+22*(diff(f1(x), x))+23*(diff(f1(x), x))+24*(diff(f2(x), x, x))+25*f2(x)+26*f2(x)+27*f3(x)+28*f3(x)+29*f3(x)+30*f3(x)+f3(x)*omega^2 = 0, diff(f1(x), x, x, x, x)+2*(diff(f1(x), x, x))+3*(diff(f2(x), x, x, x))+4*(diff(f2(x), x))+5*(diff(f3(x), x, x, x))+6*(diff(f3(x), x))+7*f1(x)+f1(x)*omega^2 = 0, 8*(diff(f2(x), x, x, x, x))+9*(diff(f2(x), x, x))+10*f2(x)+11*(diff(f1(x), x, x, x))+12*(diff(f1(x), x))+13*(diff(f3(x), x, x))+14*f3(x)+f2(x)*omega^2 = 0, f1(0) = 0, f1(1) = 0, f2(0) = 0, f2(1) = 0, f3(0) = 0, f3(1) = 0, ((D@@1)(f1))(0) = 0, ((D@@1)(f1))(1) = 0, ((D@@1)(f2))(0) = 0, ((D@@1)(f2))(1) = 0, ((D@@1)(f3))(0) = 0, ((D@@1)(f3))(1) = 0, ((D@@2)(f3))(0) = 0, ((D@@2)(f3))(1) = 0};
sys,bcs:=selectremove(has,dsys3,diff);
extra_bcs:={seq(seq(seq((D@@i)(cat(f,j))(a),j=1..2),i=2..3),a=0..1),seq(seq((D@@i)(f3)(a),i=3..5),a=0..1)};
nops(extra_bcs);
for b in extra_bcs do
  try
    print(b=-1);
    res[b]:=dsolve(dsys3 union {b=-1}, numeric,initmesh=1024,approxsoln=[omega=1,f1(x)=x*(1-x),f2(x)=sin(Pi*x),f3(x)=x*(1-x)]);
  catch:
    print(lasterror);
  end try
end do;

eval(res);
RES:=res[(D@@4)(f3)(0)];
RES(0);
plots:-odeplot(RES,[[x,1000*f1(x)],[x,1000*f2(x)],[x,f3(x)]],0..1);
#### From RES(0) you get that omega = 948.876522424828
####Notice that f1 and f2 have been multiplied by 1000 in the plot.

@Preben Alsholm 
Continuing with your new dsys3 this works for one of the extra bcs:

sys,bcs:=selectremove(has,dsys3,diff);
extra_bcs:={seq(seq(seq((D@@i)(cat(f,j))(a),j=1..2),i=2..3),a=0..1),seq(seq((D@@i)(f3)(a),i=3..5),a=0..1)};
nops(extra_bcs);
for b in extra_bcs do
  try
    print(b=-1);
    res[b]:=dsolve(dsys3 union {b=-1}, numeric,initmesh=1024,approxsoln=[omega=1,f1(x)=x*(1-x),f2(x)=sin(Pi*x),f3(x)=x*(1-x)]);
  catch:
    print(lasterror);
  end try
end do;

eval(res);
RES:=res[(D@@4)(f3)(0)];
RES(0);
plots:-odeplot(RES,[[x,1000*f1(x)],[x,1000*f2(x)],[x,f3(x)]],0..1);

@9009134 dsolve/numeric will allow a parameter in the system if provided with an extra condition that will determine that parameter.
I think the idea is that the extra parameter (omega, say) internally is represented by the ode diff(omega(t),t) =0.
The extra condition will of course be put on the unknown functions (in your case f1,f2,f3) not on omega.

I tried the following on your latest dsys3:
extra_bcs:={seq(seq(seq((D@@i)(cat(f,j))(a),j=1..2),i=2..3),a=0..1),seq(seq((D@@i)(f3)(a),i=3..5),a=0..1)};
nops(extra_bcs);
###
for b in extra_bcs do
  try
    print(b=1);
    res[b]:=dsolve(dsys3 union {b=1}, numeric,initmesh=1024);
  catch:
    print(lasterror);
  end try
end do;

As you will see I didn't have success for any of the 14 choices of extra conditions. (Originally I claimed erroneously that there were 8: there are 14 of the kind I talked about). "matrix is singular" was the error message in each case.

@erik10 I upgraded my Windows 8.1 to Windows 10 a couple of days ago. Maple 2015.1 using worksheet mode and 1-D input works without any problems according to my 2 days of experience.

@9009134 This is a totally different system!

You may find inspiration in this simple example, which you can easily solve symbolically, but here it is done numerically:

restart;
ode:=diff(y(x),x,x)+lambda*y(x)=0;
res:=dsolve({ode,y(0)=0,y(1)=0,D(y)(0)=1},numeric);
res(0);
evalf(Pi^2);
plots:-odeplot(res,[x,y(x)],0..1);
res2:=dsolve({ode,y(0)=0,y(1)=0,D(y)(0)=1},numeric,approxsoln=[lambda=40,y(x)=sin(2*Pi*x)]);
res2(0);
plots:-odeplot(res2,[x,y(x)],0..1);

Your problem is to find a nonzero boundary condition that works for your case. You have plenty to try. If you choose the extra condition to be of the form (D@@k)(f_i)(a) = 0 where k is less than the order of f_i and where a is 0 or 1, I think you have 2*(4+4+6) - 14 = 14 different ones. (I corrected an error in this calculation)!!

@baustamm1 Yes I meant building blocks in that sense.

Whar are the intial and boundary conditions here in this case?

@Carl Love Your original version works in Maple 2015.1 with 1-D math input.

@baustamm1

The output from pdsolve on your system sys1 is not to be considered the general solution, but rather as building blocks for any solution as when using separation of variables.

You could try doing
infolevel[pdsolve]:=3;
before executing the pdsolve command.

@Axel Vogt No I don't think so.
Try setting
infolevel[dsolve]:=3:
Then when executing the dsolve statement you will see the procedure for evaluating YP[1] is
proc (N, X, Y, YP) option `[Y[1] = f(x)]`; YP[1] := ff(X); 0 end proc;
Thus ff(X) remains unevaluated. Thus no problem.
The problem appear right away when defining the ode if you didn't have the return unevaluated part, as in
ff1:= x -> evalf(Int(exp(t), t= 0 .. x, method = _d01ajc));
ode1:=D(f)(x)=ff1(x);
You will notice the disappearance of option and evalf before we even get to the task of solving.

@Mac Dude You should be able to use Physics:-Assume without loading the package. Just use the long name Physics:-Assume.
Example:
restart;
Physics:-Assume(a>0);
sqrt(a^2);
about(a);
Physics:-Assume(a=a);
about(a);

@Kitonum If the order is important then notice that your solution doesn't preserve the order, even in the example given.
Here is another example:

restart;
u:=[x,z,y]; v:=[b,a,x];
[op({op(u), op(v)})];
ListTools[MakeUnique]([u[],v[]]);

OK, you say that if you want the result sorted then etc., but sorted according to which criterion? In recent versions we get lexicograhic ordering when using sets, thus 'a' preceeds 'b' etc.

@acer I just submitted an SCR with title "Bug in indets".

Maybe superfluous to mention, but you have to remember to set Typesetting to Extended in the options menu.

The problem is the same for int and Int. The remedy by acer still works.

restart;
w:=[int(f(x),x),Int(f(x),x)];
indets(w,specfunc(int));
indets(w,specfunc(Int));
convert(w,Int);
convert(w,int);
value(w);
indets(w,specfunc(int) &under (convert,int)); #Empty
indets(w,specfunc(int) &under value);#Empty
indets(w,specfunc(Int) &under (convert,Int));#Empty
## acer's remedy of converting global as the last action:
indets(w,specfunc(int) &under (convert,compose,int,`global`));
indets(w,specfunc(int) &under (rcurry(convert,`global`)@value));
indets(w,specfunc(Int) &under (convert,compose,Int,`global`));
## The converting global is essential:
indets(w,specfunc(XX) &under (curry(subs,int=XX)@value));
indets(w,specfunc(XX) &under (rcurry(convert,`global`)@curry(subs,int=XX)@value));
###########################################################
The problem is indeed rather confusing. Yet another example, this time with f and %f (and in parallel f(x) and %f(x)):
restart;
w:=[f,%f];
w(x);
value(%f);
indets(w,identical(f));
indets(w,identical(%f));
value(w);
indets(w,identical(f) &under value);#Empty
indets(w(x),specfunc(f) &under value);#Empty
indets(w,identical(%f) &under curry(subs,f=%f)); #OK
indets(w(x),specfunc(%f) &under curry(subs,f=%f)); #OK
indets(w,identical(f) &under (rcurry(convert,`global`)@value)); #Surprise!
indets(w(x),specfunc(f) &under (rcurry(convert,`global`)@value)); #OK!
indets(w,identical(%f) &under (rcurry(convert,`global`)@curry(subs,f=%f))); #Surprise!
indets(w(x),specfunc(%f) &under (rcurry(convert,`global`)@curry(subs,f=%f))); #OK!