for k = 1:101
    x(k) = (k-1)*0.01;
    y(k) = (k-1)*0.01;
end
for p = 1:101
    for q = 1:101
        u00(p,q) = 0;
        u01(p,q) = 0;
        u02(p,q) = 0;
        u03(p,q) = 0;
 %
 % NOTE: In this simple example the solution contains a single mode
 % of n = 2 and m = 3, therefore the summation in the following
 % do-loop is only from n = 2-2 and m = 3-3.  In general, the
 % summation should be from n = 1-N, m = 1-M, where N and M are
 % the wavenumbers in x- and y-direction at which we truncate
 % our infinite series
 %
        for n = 2:2
        for m = 3:3
 %
 % In a more complicated problem amn will not be just 1 but
 % will depend on m and n;  Modify the following line accordingly
 %
            amn = 1;
 % solution for t = 0
            t = 0;
            ct = cos(sqrt(m*m+n*n)*pi*t);
            u00(p,q) = u00(p,q)+amn*sin(n*pi*x(q))*sin(m*pi*y(p))*ct;
 % solution for t = 45/(100*sqrt(13))
            t = 45/(100*sqrt(13));
            ct = cos(sqrt(m*m+n*n)*pi*t);
            u01(p,q) = u01(p,q)+amn*sin(n*pi*x(q))*sin(m*pi*y(p))*ct;
 % solution for t = 55/(100*sqrt(13))
            t = 55/(100*sqrt(13));
            ct = cos(sqrt(m*m+n*n)*pi*t);
            u02(p,q) = u02(p,q)+amn*sin(n*pi*x(q))*sin(m*pi*y(p))*ct;
 % solution for t = 100/(100*sqrt(13))
            t = 100/(100*sqrt(13));
            ct = cos(sqrt(m*m+n*n)*pi*t);
            u03(p,q) = u03(p,q)+amn*sin(n*pi*x(q))*sin(m*pi*y(p))*ct;
        end
        end
    end
end
subplot(2,2,1);
pcolor(x,y,u00)
shading interp
hold on
for m = 1:9
    pm = m*0.1;
    contour(x,y,u00,pm,'k-')
    nm = m*(-0.1);
    contour(x,y,u00,nm,'k--')
end
hold off
subplot(2,2,2);
pcolor(x,y,u01)
shading interp
hold on
for m = 1:9
    pm = m*0.1;
    contour(x,y,u01,pm,'k-')
    nm = m*(-0.1);
    contour(x,y,u01,nm,'k--')
end
hold off
subplot(2,2,3);
pcolor(x,y,u02)
shading interp
hold on
for m = 1:9
    pm = m*0.1;
    contour(x,y,u02,pm,'k-')
    nm = m*(-0.1);
    contour(x,y,u02,nm,'k--')
end
hold off
subplot(2,2,4);
pcolor(x,y,u03)
shading interp
hold on
for m = 1:9
    pm = m*0.1;
    contour(x,y,u03,pm,'k-')
    nm = m*(-0.1);
    contour(x,y,u03,nm,'k--')
end
hold off