@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.

@J4James I just tried to google  elliptic pdes in Maple  and one item which came up looked promising:

http://www.math.unb.ca/~seahra/4503/maple_2011.html
The author is Sanjeev Seahra.
You can get a worksheet and a pdf file from that page on Linear Elliptic PDEs.
I haven't taken the time to look at it, but just seen that it is available and can be downloaded without any problem.

@wo0olf You can use sol:-value with the option output=listprocedure like this:

Tv:=subs(sol:-value(output=listprocedure),T(x,y));
Tv(.01,.00002); #Example of use
100/Tv(.01,.00002); #h at the same (x,y) point

Your system is not an autonomous system, i.e. the right hand sides depend on t explicitly through the function sin(2*Pi*t). (I'm assuming that by sin2Πt you mean sin(2*Pi*t)). Thus the vector field associated with the odes for S(t) and i(t) is time dependent (I use 'i' instead of 'I', because I is the imaginary unit in Maple).
Therefore you ought to tell us in more detail what kind of picture you would like to see.

@Chia No, as an example, try:
pde:=diff(u(x,y),x,x)+diff(u(x,y),y,y)=0;
bcs:={u(x,0)=0,u(x,1)=0,u(0,y)=y,u(1,y)=y^2};
pdsolve(pde,bcs,numeric);
The error message is:

Error, (in pdsolve/numeric) unable to handle elliptic PDEs

In the help page for pdsolve/numeric the solver is described as time based, which I'm sure means that the algorithm "marches" forward in time (t) from known initial conditions all along satisfying the boundary conditions at each end of the x-interval. You get no complaints therefore if you impose initial and boundary conditions which allow that to happen:
ibcs:={u(x,0)=0,D[2](u)(x,0)=0,u(0,y)=y,u(1,y)=y^2};
res:=pdsolve(pde,ibcs,numeric,timestep=.001,spacestep=.001);
Whether it is a reasonable problem (ill-posed, or interesting) is another matter.

This problem was discussed in an answer to your previous question:

http://www.mapleprimes.com/questions/202074-Pde-Solver-Problem

Just to give an answer:

res:=u(x, t) = Sum((Int(f(x)*cos(Pi/4*(2*n+1)*x), x = 0 .. 2))*cos(Pi/4*(2*n+1)*x)*exp(-(Pi/4*(2*n+1))^2*a^2*t), n = 0 .. infinity);

This is found by hand the same way Maple successfully solves the problem with the boundary conditions
u(0,t) = 0, u(2,t) = 0.