@farahnaz You have changed the system considerably, I notice.

Starting with the new dsys3, which I shall not print here) I did:

sys, bcs := selectremove(has, dsys3, diff);
indets(sys,specfunc(diff)); #Finding the orders
solve(sys,{diff(f1(x),x$3),diff(f2(x),x$4),diff(f3(x),x$4),diff(f4(x),x$2)}); #Try isolation of highest derivatives: FAILURE
newsys:={sys[1],diff(sys[2],x),sys[3],sys[4]}: #sys[2] is differentiated
indets(newsys,specfunc(diff)); #Same as for sys
solve(newsys,{diff(f1(x),x$3),diff(f2(x),x$4),diff(f3(x),x$4),diff(f4(x),x$2)}):
nops(%); #answer 4, thus highest derivatives can be solved for
##Because of the differentiation of sys[2] a constraint given by sys[2]=0 is to be kept.
## It is enough to require it satisfied at x=1 (or x=0) since sys[1] must be constant
bcs2:=eval[recurse](convert(sys[2],D),{x=1} union bcs);
nops(bcs union {bcs2}); #Answer 13, i.e. the same as the sum of the highest orders 3+4+4+2.
##However, there is the parameter omega. So we need an extra condition.
##Simple possible ones to examine are
extra_bcs := {seq(seq(seq(((D@@i)(cat(f, j)))(a), j = 2 .. 3), i = 2 .. 3), a = 0 .. 1), seq(seq(((D@@i)(f1))(a), i = 1 .. 2), a = 0 .. 1), seq(seq(((D@@i)(f4))(a), i = 1 .. 1), a = 0 .. 1)};
##You should definitely replace omega^2 by omega2 (i.e. a name):
newsys2:=subs(omega^2=omega2,newsys):
## I had immediate luck with just one of these (notice that abserr=2e-6):
for b in extra_bcs do
  try
    print(b=1);
    res[b]:=dsolve(newsys2 union bcs union {bcs2} union {b=1},numeric,method=bvp[middefer],initmesh = 128, approxsoln = [omega2 = 1, f1(x) = x*(1-x), f2(x) = sin(Pi*x), f3(x) = x*(1-x), f4(x) = x^2*(1-x)^2],abserr=2e-6)
  catch:
    print(lasterror)
  end try
end do;
#LUCK with
RES:=res[(D@@3)(f2)(1)];
##To find omega^2 = omega do
RES(0);
plots:-odeplot(RES,[x,f1(x)],0..1); #Example plot

I tried setting Digits:=16.  Computing sol(300) used 44.5 sec and my matrix was computed in 474 sec. Thus the ratio is more favorable for me, but still rather large at roughly 11.
My procedure uses too many symbolic constructs to make significant use of evalhf.
I started writing it a few years ago for my own educational purpose. It has no error control, so my use of the adjective "amateurish" was not just a show of false modesty.

@Maliha Saleem In the worksheet downloaded to my computer the epsilon in the pde 'ge' is 'varepsilon', which prints the same as epsilon, but you assign to epsilon later, not to varepsilon.
This has nothing to do with the errors at the bottom, which (as tomleslie has pointed out) can be fixed by converting to 1D-math input.

@Maliha Saleem I work exclusively in Worksheet mode and 1D -math input.

You are apparently using 2D input and the pop up dialogue asks whether you want function definition or remember table assignment. So you are the one to decide.
In 1D input you are not faced with that choice: If you do u[0](zeta):=whatever, you get a remember table assignment.
If I want a function definition then in this case I would do:
u[0]:=unapply(solve(eq[0],u[0](zeta)),zeta);
You can do that too in 2D.

More interesting: Why do you want to assign to u[0] or u[0](zeta) in the first place?

With the code as you added it later I still don't have any problems:

restart;
bc1 := u(x, y)-0;
dydxf := (1/2)*(-u[m+2](zeta)-3*u[m](zeta)+4*u[m+1](zeta))/h;
bc1 := subs(diff(u(x, y), x) = subs(m = 0, dydxf), u(x, y) = u[0](zeta), x = 0, bc1);
eq[0] := bc1;
u[0](zeta):=solve(eq[0],u[0](zeta));

In a fresh worksheet I had no problem with this:

u[0](zeta):=0;

so obviously we need to know more. Give us the whole code starting with restart or upload a worksheet.

What impresses me with example 2 is the speed. I compared with my homemade and thus amateurish DDE solver. I'm happy that it produced the same graphs as dsolve did, but as I expected it took quite a lot longer: dsolve on my computer took roughly 5 seconds to compute sol(300), whereas my program which produced a 30201x7 matrix took approximately 70 times as long.

@Rouben Rostamian  Yes, you can use an approximate solution. In this case we know the solution so the following is cheating to some degree:
##Verifying an exact solution
odetest(y(x)=1/2*sin(2*x),eval({de,bc},omega2=4));
## I'm not using the exact solution but something vaguely similar, and it happens to work:
res:=dsolve({de,bc}, numeric,approxsoln=[y(x)=sin(2.4*x),omega2=3]);

@kyanashayan If possible include the integral in its differentiated form in the system.
I will give a simple example where I use two different methods:

restart;
ode:=diff(y(x), x, x) = -(8*omega*(1-exp(-8*x))*exp(-8*x)/(2*x^2)-Pi^2/(32*x))*y(x);
bcs:=y(0)=0,D(y)(0)=1,y(1/2)=0;
## First method using 'output=listprocedure':
res:=dsolve({ode,bcs},numeric,method=bvp[middefer],output=listprocedure);
plots:-odeplot(res,[x,y(x)],0..1/2);
Y:=subs(res,y(x)); #Fishing out the procedure for computing y(x)
A:=Int(sin(x*y(x))+cos(x),x=0..1/2); #The simple example of an integral
evalf(eval(A,y=Y)); #Replacing y by Y and using numerical integration
## Second method for the same integral:
ode2:=diff(a(x),x)=sin(x*y(x))+cos(x);
res2:=dsolve({ode,ode2,bcs,a(0)=0},numeric,method=bvp[middefer]);
res2(0.5);

@Rouben Rostamian  Replacing omega^2 by omega2 solves the problem.

@farzane There is no way that you can get symbolic results in that case. You will have a characteristic polynomial of degree 8 containing 30 (!!!) parameters. As there is no general formula for the roots of a polynomial of degree 5 and above you stand absolutely no chance of getting a symbolic result.

Incidentally the exponential function exp(x) is written like that, NOT like this e^x.
You can correect that after the fact by doing as follows:

indets(Sys,`^`); #Just to see the problem
op(e^x); #The operands of e^x
Sys2:=evalindets(Sys,`^`,s->exp(op(2,s)));
indets(%,function);
Var:=[b1,c1,a2,b2,c2,d2,a3,d3];
indets(Sys2,name); #Notice the indexed v
nops(%) - 8; ## 30 (!!!)
(A, w) := LinearAlgebra[GenerateMatrix]( Sys2, Var ); #w :You cannot use v as it appears in A

op(1,A); #8x8

@al-faraj Did you try the code I gave? If you did you should have received the answer

10.11073088

And I don't think that I misunderstood your code.
I wasn't aware of a Danish version of Maple.
In your Maple version is the decimal separator (or whatever the term is) a comma?
What comes out of evalf(Pi); in your version
3.141592654
or
3,141592654

?

@maple fan The branches for a^z are those which are consequences of its definition:

a^z = exp(z*log(a)), where log is any branch of the logarithm.

Those are

log(z) = ln(abs(z)) + I*n*argument(z)

where n is any integer and where argument(z) is any argument for z.

@mskalsi Although the use of the files maple.hdb, maple.ind, maple.lib some (many) versions ago was replaced by archive files ending in .mla the use of the old method still works.
I cannot test it on the ones you mention since it seems you have to pay to get access to them.
However, I happen to have some old maple.ind, maple.lib (no maple.hdb). They contain among other things a package called DMat.

On a computer with Maple 2015 I did
libname:="E:/DMat/libDMat",libname;
with(DMat);
##The content showed up.
## In libDMat are just the files maple.ind and maple.lib

After that on another computer with Maple 12 I did the similar thing (the path is different). It worket too.

So I cannot tell you why you cannot make it work.

#Just an afterthought: Is Rif from the link really a package or is it just a collection of procedures not made into a package. If that is so then with(Rif); would result in an error.
Another thing: Is the "R" in Rif really capitalized? Should it be rif?
Thirdly, as Markiyan is mentioning, rif is implemented in Maple. In Maple 12 (I don't have access to Maple 13) when I do
?rif
up comes a help page "Overview of the Rif Command Set Version 1.1"  (just as in Markiyan's comment).

@mskalsi Are you sure that the path is correct?
My path
"E:/RIFfolder"
was just an example. You need to give the exact path to the place where you put the Rif-files.
(In particular: It doesn't seem likely that the drive is E as in my example).