@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

To get a useful response you will have to provide the details. You could upload a worksheet or give us the full text of the relevant code.

You got an integral involving the unknown w(x,t) in PDE1. So the system PDE1, PDE2 does not only include derivatives but also an integral, so it is certainly not obvious to me what you could do.

@emrantohidi No ready made code is available in Maple 18 or earlier for the numerical handling of your equations.

@emrantohidi It works fine in Maple 18. You should replace the constant e in the boundary condition
u(x, 0, z, t) = e*exp(3*t+x+z)
with
u(x, 0, z, t) = exp(1)*exp(3*t+x+z)
or (of course)
u(x, 0, z, t) = exp(3*t+x+z+1).
e in Maple is just another name.

I just checked the code in Maple 15 also. It came up with the solution. The code:

restart;
PDE1 := diff(u(x, y, z, t), t)-(diff(u(x, y, z, t), `$`(x, 2)))-(diff(u(x, y, z, t), `$`(y, 2)))-(diff(u(x, y, z, t), `$`(z, 2))) = 0;
IC := {u(x, y, z, 0) = exp(x+y+z)}; BC := {u(0, y, z, t) = exp(3*t+y+z), u(1, y, z, t) = exp(3*t+y+z+1), u(x, 0, z, t) = exp(3*t+x+z+1), u(x, 1, z, t) = exp(3*t+x+z+1), u(x, y, 0, t) = exp(3*t+x+y), u(x, y, 1, t) = exp(3*t+y+x+1)};
Sol := pdsolve({PDE1} union IC union BC);

However, in Maple 12 neither the 2-dim nor the 3-dim works.
I don't have Maple 13 or 14 available here.

I use a Danish keyboard, but admittedly I always use the worksheet interface and 1D (aka Maple) input, so my reply may not be relevant.
However, I just tried the keystrokes [CTRL]+[ALT]+0 in Maple 17 and 18 and it produced the }  just as [ALT-GR]+0 did.
When I try  [CTRL]+= the input just gets very small.

@emrantohidi Your link refers to a long paper and your question seems to be of a general nature.
I'm not sure that I would be the right person to answer that. In any case I would certainly suggest that you ask the question as a new question and not here buried in some other question. That should give it a chance of being seen by more people.
Also I strongly suggest that you provide a specific example of the type of problem you are talking about.
The more concrete the better.