@Joe Riel I only got Iterator.hdb and Iterator.mla in Iterator/lib.
The error message from with(Iterator);  is

Error, invalid input: with expects its 1st argument, pname, to be of type {`module`, package}, but received Iterator

In Maple 17, where the help works I get the same error message.

The version I installed (or thought I installed) is 2.3.4.

@Joe Riel I got the most recent version of the Iterator package from the Maple Application Center.
I have used an earlier version before. I'm using Windows Vista. In Maple 18 the libname seems quite right:
"C:\Program Files\Maple 18\lib", "C:\Users\alsholm\maple\toolbox\Iterator\lib", "."
That User directory is the place where Iterator.hbd and Iterator.mla are located.
However, ?Iterator doesn't bring up the help for Iterator, but something from combstruct. Also with(Iterator); just gives an error message.
In Maple 17, however, where libname still includes the same User directory, ?Iterator does bring up the help page for Iterator. But also here with(Iterator); doesn't work.
Any guess as to what I might be doing wrong?

@nm It is unsafe to use sum since in that case Maple will only execute the procedures a and b twice: once for n=0 and once for general n. The exceptional cases will not be caught by int. If you replace sum with add the procedures a and b will be called maxN + 1 times so that int will work on concrete values of n.
To see this you can insert a print statement in the procedure a:
a:=proc(n) print(n);
           int(f(x)*cos(n*freq*x),x=from_..to_) /denomC;
         end proc;

Your equation is missing two multiplication signs, they have been put in by DJJerome below.

@Chaimongkol I entered:
ode:=(4*Wi^2*t^3+4*t+4*sqrt(1-t^2)*De)*u(t)+(-8*sqrt(1-t^2)*De*t-4*L-4*De^2*t^3)*diff(u(t),t,t)=0;
infolevel[dsolve]:=5:
dsolve(ode);
And got a lot of printout. At one point it said:
-> Trying a Liouvillian solution using Kovacic's algorithm
   <- No Liouvillian solutions exists

Then it continued with something else, but didn't finish with anything before my patience ran out.

@baharm31 Saving w-values to a vector or matrix is easily done. Here is the simple example I used above.

restart;
pde:= diff(w(x, t),x,x)-diff(w(x,t),t)=0;
pi:=evalf(Pi): #Instead of just 3.14.
res:=pdsolve(pde,{w(x,0)=sin(x),w(0,t)=0,w(pi,t)=0},numeric);
res:-plot3d(t=0..1);
W:=subs(res:-value(output=listprocedure),w(x,t));
W(pi,0.1); #test
#To see what is going on m and n are fairly small.
#Increasing x- and t-values along rows and columns respectively.
t0:=0: t1:=1: a:=0: b:=pi: m:=10: n:=8:
h:=(b-a)/(n-1): k:=(t1-t0)/(m-1):
A:=Matrix(m,n,(i,j)->W(a+h*(j-1),t0+k*(i-1)));
A[..,3]; #The third column: fixed x
A[4,..]; #The fourth row: fixed t
plots:-matrixplot(A);

@baharm31 I tried to come up with an example with an integral where an exact solution to the problem is known.
Here is one which can be varied to some degree.
I determine ft such that u(x,t) = sin(x+t) is a solution to the pde and the initial and boundary conditions.
restart;
pde:=diff(u(x,t),t,t)=diff(u(x,t),x,x)+int(diff(u(x,t),x)*u(x,t),x=0..1/2)+ft;
eval(pde,u(x,t)=sin(x+t));
S:=solve(%,{ft});
pde:=subs(S,pde);
ibcs:={u(x,0)=sin(x),D[2](u)(x,0)=cos(x),u(0,t)=sin(t),u(1,t)=sin(1+t)};
res:=pdsolve(pde,ibcs,numeric,timestep=0.05,spacestep=0.05);
res:-animate(u(x,t)-sin(x+t),t=5);
res:-plot3d(t=5);
plot3d(sin(x+t),x=0..1,t=0..5);
res:-plot3d(u(x,t)-sin(x+t),t=5);
#The maximal discrepancy is about 4 percent as you can see.
Variations: As integrand you can try diff(u(x,t),x)^2, u(x,t)^2, u(x,t)^3, which work, but it is puzzling that u(x,t) gives an error:
Error, (in pdsolve/numeric/process_PDEs) selecting function must return true or false

The same error message comes with the integrand diff(u(x,t),x).
################################ Comment #####################
I tried a similar example using my own unsophisticated code which handles pdes of the form
diff(u(x,t),t) = f(x,t,u,ux,uxx)
where f is just a function of the 5 real variables (not functions) x,t,u,ux,uxx. The unamended code delivered a result on an f that included an integral of the function u(x,t)^2 w.r.t. x which was not at all ridiculous.
I conclude that the fact that pdsolve occasionally seems to handle integrals is unintended. Also there is no claim made in the help pages that it should do that.
The integral int(u(x,t)^2,x=0..c) is in my unamended code just replaced by c*u[j,n]^2. That clearly is not a very good substitute for the integral. It ought to be something like h*add(u[jj,n]^2,jj=0..j1), where j1 corresponds to the x-value c and where h is the spacestep.

@baharm31 I tried setting infolevel[int]:=5.
It produced a lot of printout, but ended claiming to have succeeded (to my great surprise).
Then I tried replacing the integral with 0 and there clearly v was just identically zero, which it ought to be.
After that I tried replacing the integral by its inegrand:
pde20:=evalindets(pde2,specfunc(anything,int),s->op(1,s));
res2:=pdsolve({pde1,pde20,pde3},bcs,numeric);
res2:-plot3d(w(x,t),t=0..tmax);
res2:-plot3d(v(x,t),t=0..tmax);
res2:-plot3d(u(x,t),t=0..tmax);
p:=res:-value(v);
p(xmax*.123456,tmax*.8765);
p2:=res2:-value(v);
p2(xmax*.123456,tmax*.8765);
#Clearly different, so maybe the integral is handled OK, but I'm still sceptical.
This time I haven't tried Dirac, but from what I said earlier that would have to be treated like 0.

@baharm31 After debugging pdsolve/numeric it appears that a check in the procedure `pdsolve/numeric/par_hyp` of the orders of the highest derivatives w.r.t. the space variable is performed.
It seems that it is necessary for the algorithm that the sum of the orders for the dependent variables (v and w) should be the same as the sum of the highest orders for each equation. In your case the first sum is 0+4 = 4 and the second sum is 2+4 = 6. Since 4<>6 you get the error message about not being close to standard form.
To make a small test of this try replacing diff(w(x,t),x,x) by w(x,t) in PDE1a. This will bring down the second sum to 0+4= 4, which makes it agree with the first:
res:=pdsolve({subs(diff(w(x,t),x,x)=w(x,t),PDE1a),PDE2a},bcs2,numeric);
res:-plot3d(w(x,t),t=0..tmax);
res:-plot3d(v(x,t),t=0..tmax);

#############Introducing a third dependent variable ###############
It appears that by introducing one more dependent variable the algorithm will work:
pde1:=diff(w(x,t),x,x)=u(x,t);
pde2:=subs(pde1,PDE1a);
pde3:=subs(pde1,PDE2a);
#Notice that the first sum is now (for u,v,w): 2+0+2 = 4.
#The second sum is (pde1,pde2,pde3): 2+0+2 = 4.
bcs:={v(x, 0) = 0, w(0, t) = 0, w(x, 0) = 0, D[1](w)(0, t) = 0, D[2](w)(x, 0) = 0, u(.57, t) = 0, D[1](u)(.57, t) = 0};
res:=pdsolve({pde1,pde2,pde3},bcs,numeric);
res:-plot3d(w(x,t),t=0..tmax);
res:-plot3d(v(x,t),t=0..tmax);
res:-plot3d(u(x,t),t=0..tmax);

@baharm31 bcs2 has the initial and boundary conditions for both of the pdes. Thus the syntax ought to be:

PDES := pdsolve({PDE1a,PDE2a}, bcs2,numeric);

However, Maple gives the error:

Error, (in pdsolve/numeric/par_hyp) input system is too far from a 'standard' form (see ?pdsolve,numeric for more detail)

Then when you go to that help page it says in the section with the heading 'Time-based Solver':

"The first mode of operation uses the default method, which is a centered implicit scheme, and is capable of finding solutions for higher order PDE or PDE systems. The PDE systems must be sufficiently close to a standard form for the method to find the numerical solution."

That doesn't really give us the promised "more detail".

@baharm31 You will have to get rid of the integral, but then you have an entirely different system.
To solve these two much simplified equations together numerically you have to replace v(t) by v(x,t) in both equations, which is easily done by subs:
PDE1a:=subs(v(t)=v(x,t),PDE1); #Similarly in PDE2
Secondly, you have to provide an initial condition for v(x,t) of the form v(x,0)=f(x), where f is known.
Then simply do
res:=pdsolve({PDE1a,PDE2a], bcs union {v(x,0)=f(x)},numeric); #Add options as you prefer

@baharm31 By without Dirac I thought that you meant to leave out the term having Diract in it. That term is the only one involving v.
So what do you mean by PDE2 without Dirac?

@baharm31 You said that I could ignore the Dirac part. That makes PDE2 a pde in w(x,t) only. Thus that one can be solved numerically independent of PDE1.
Now PDE1 involves the integral of derivatives of w. Thus you need to get those from the solution of PDE2 before you can solve PDE1. That was what I was getting at with a somewhat simpler example: How to get the derivatives of w that we need.

@Chia You applied fsolve to Re(t). BesselJ(0,y) is real for real values of y. Thus the real part of t is
evalc(Re(t)) assuming y>0;
i.e. the real part is also zero where cos(y) = 0. Thus you should actually get 6 roots for Re(t) for y in 0..10:
Illustration:
plot([cos(y),BesselJ(0,y)],y=0..10);
plot(Re(t),y=0..10);
So, yes fsolve is not perfect.

Now RootFinding:-Analytic you applied to t on a rectangle of complex y-values. It found 7 roots, but only 3 of these were real: all 3 were zeros of BesselJ(0,y) as also found by BesselJZeros(0.,1..3);

@baharm31 It seems that you cannot get values for derivatives like diff(w(x,t),x,x) or diff(w(x,t),t) from the use of value, although when one looks at the response from trying to do so it appears possible.
So the situation is different from the corresponding one for odes.
Notice though that there is no claim made in the help page for pdsolve/numeric that it should be possible.
Consider for clarity this simpler example:

restart;
pde:= diff(w(x, t),x,x)-diff(w(x,t),t)=0;
res:=pdsolve(pde,{w(x,0)=sin(x),w(0,t)=0,w(3.14,t)=0},numeric);
res:-plot3d(t=0..1);
res:-plot3d(diff(w(x,t),x),t=0..1); #Doesn't work: flat zero
res:-value(diff(w(x,t),x),output=listprocedure); # Appears to work
%(2,.5); #But doesn't work: 0
#So we try using fdiff on this procedure:
W:=subs(res:-value(output=listprocedure),w(x,t));
W(2,.3);
Wx:=proc(x,t) fdiff(W,[1],[x,t]) end proc; #Computes wx for given values of x and t.
Wx(2,.3);
plot3d(Wx,0..3,0..1);
#Try computing an integral
p:=proc(t) Wx(2.345,t)*W(2.345,t) end proc;
p(0.5);
evalf(Int(p,0..1,epsilon=1e-9)); #Work OK