I don't know what happens here. As I told you I have an old version of Mathematica, which doesn't understand many of your commands.
Here is a code in which I changed the length of the square to 3.6. As before the flow of heat at the boundaries is zero (that are your bc's) and the q gives a constant supply of energy at each point, giving rise to a uniform rise of T throughout the square (meaning. you don't need your partial diffequation)
ec = rho*cp*D[T[x, y, t], t] == k*(D[T[x, y, t], x, x] + D[T[x, y, t], y, y]) + q;
f = 6780000; e = 30000; \[Epsilon] = 55.33*10^-12;
k = 104; h = 10; rho = 410; cp = 2.2; To = 300; Tamb = 323;
q = f*e^2*\[Epsilon];
alpha = k/(rho*cp);
ci = T[x, y, 0] == To;
(* new Length of sides of square*)
xL = 3.6;
cc = {D[T[x, y, t], x] == 0 /. x -> 0, D[T[x, y, t], y] == 0 /. y -> 0,
D[T[x, y, t], x] == 0 /. x -> xL, D[T[x, y, t], y] == 0 /. y -> xL};
sln = NDSolve[{ec, ci, cc}, T, {t, 0, 20}, {x, 0, xL}, {y, 0, xL}];
fT = T /. Flatten[sln]
Animate[Plot3D[fT[x, y, t], {x, 0, xL}, {y, 0, xL},
PlotRange -> {200, 9000}], {t, 0, 20}]
My email is h.dolhaine@gmx.de