@Preben Alsholm Thanks a lot.

Now, the code runs and It is able to calculate and plot.

However, the calculation in Maple seems not correct!  For example, I gave center points as (x,y)=(0.11143,0.94020) in Gaussian function which must have a value 5 (Amplitute) i.e r(r,phi)=f(0.9415,1.4528)=5  (polar coordinates (r,phi)=(0.9415,1.4528)) and computed w(0.9415,1.4528,0)=3.561804477 i.e the solution at time=0. This must be the exactly same value with the gaussian function (initial value) f(0.9415,1.4528).

This might be because of the huge peak value in Gaussian function that numerical integration cannot handle. 

If I decrease sigma(x and y spreads of the blob) i.e no sharp peak, then f(r',phi')=w(r',phi',0).

Any idea(s) what kind of methods in Maple should be used for such huge gradients?

Download replied_2d_soln.mw

@mehdi jafari Thanks a lot. I got the idea. I tried the following

f(x,y):=x+y:
c:=10^(-1): G(x,y,xi,eta,t):=(1+2*(sum((exp(-Pi^(2)*n^(2)*c*t))*cos(n*Pi*x)*cos(n*Pi*xi),n=1..10)))*(1+2*(sum((exp(-Pi^(2)*m^(2)*c*t))*cos(m*Pi*y)*cos(m*Pi*eta),m=1..10))):
w(x,y,t):=int(int(f(xi,eta)*G(x,y,xi,eta,t), eta=0..1), xi=0..1):
plot(w(.51086, .49427, t), t = 0 .. 1);
for t from 0 to 100 do w(.51086, .49427, 0.01*t) end do;

by doing that, i can get 101*2 matrix.

Thanks again.

I used for loop

for t from 0 to 100 do w(.51086, .49427, t), t*0.01) do

to get the y_coordinates. Hence I can write 101*2 Martix. 

@mehdi jafari thanks for your prompt comment. Here is the maple code;

f(x,y):=x+y:
c:=10^(-1): G(x,y,xi,eta,t):=(1+2*(sum((exp(-Pi^(2)*n^(2)*c*t))*cos(n*Pi*x)*cos(n*Pi*xi),n=1..10)))*(1+2*(sum((exp(-Pi^(2)*m^(2)*c*t))*cos(m*Pi*y)*cos(m*Pi*eta),m=1..10))):
w(x,y,t):=int(int(f(xi,eta)*G(x,y,xi,eta,t), eta=0..1), xi=0..1):
plot(w(.51086, .49427, t), t = 0 .. 1);
plottools:-getdata(%); M := %[-1]; ExportMatrix("new2.dat", M);

It gives me 200*2 matrix. For example; first 5 coordinates (x,y) is 

0.0000000000000000e+00 1.0054420984938566e+00
5.2576050251256284e-03 1.0053136811087655e+00
9.8322190452261315e-03 1.0052447682494321e+00
1.4976850954773869e-02 1.0051969485736694e+00
2.0155584824120606e-02 1.0051685044159542e+00

But I want to have the coordinates (x,y ) something like this

0.00   y0
0.01   y1
0.02   y2
...       ...
1.00   y101

For the above coordinates, I need 101*2 Matrix with spesific x coordinates.

I have realized that the command cannot be seen now. It might be better if I type here

plottools:-getdata(%); M := %[-1]; ExportMatrix("new2.dat", M)