@Christopher2222 I delete quite a few every day. It is irritating though.

@farahnaz Since you have a different system (again) I think you ought to start a separate thread for that one. That way other people will be more likely to participate. When a new problem is introduced in comments people not participating in the discussion from the start will tend to ignore it.

@Carl Love Correcting the obviously erroneous abserr in the worksheet to 1e-4 in the latter part cures the problem while keeping middefer. 
In the first part I used abserr = 3e-4 and maxmesh = 1024 instead.

You are not using an abserr of 10^(-8), which I believe you think.
If I remember correctly, I have pointed this out earlier.
The abserr you are actually using is abserr = .466507380209733.
What you are writing when you write 1.1^(-8) is in fact (1.1)^(-8), which is .466507380209733.

I don't see any problem with your worksheet, but why do you use output=Array(.....) ?
Why not leave that out entirely?
You will then have to be more careful with the range argument in odeplot. With output=Array it is ignored.

To get the full picture (without output=Array) use the range 0..10  (N=10).
To get a close-up at t=0 use e.g. 0..0.1.

@farahnaz I don't see any way that you can do this in Maple (short of coding it all yourself).

1. The boundary conditions are given on all 4 sides of a rectangle: (0,z), (x,h/2), (x,-h/2), (L,z).
Thus the requirement for pdsolve that the system be an initial and boundary value problem is not met. pdsolve uses a time based solver. Thus you can give initial conditions for time = t0 (time will be either x or z in your case) and boundary conditions on two sides perpendicular to the time axis.
2. The problem reported by the error message in your latest worksheet must be corrected: On the line t=t0 (whatever t is, x or z) you can only have terms with that t fixed (at t=t0) not the space variable. Similarly on a line with fixed space variable you can only have conditions with that fixed.
3. I see no solution to the problem I mentioned in my previous comment and exemplified.

Rest assured that even if you correct 2 out of the 3 problems above the third will prevent pdsolve from working.
pdsolve just reports the error it happens to find first and then it stops. 

@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);