Solving a 2D heat equation on a square with Dirichlet boundary conditions

ClearAll[x, y, t, f];

f = 10*Sin[Pi*x]^10*Sin[Pi*y]^10*Sin[Pi*t]^10;
pde = D[u[x, y, t], t] == D[u[x, y, t], {x, 2}] + D[u[x, y, t], {y, 2}] + f;
ic = u[x, y, 0] == 0;
bc = {u[0, y, t] == 0, u[1, y, t] == 0, u[x, 0, t] == 0, u[x, 1, t] == 0};

sol = NDSolve[{pde, ic, bc}, u, {x, 0, 1}, {y, 0, 1}, {t, 0, 4}];

Manipulate[
 Plot3D[Evaluate[u[x, y, t0] /. sol], {x, 0, 1}, {y, 0, 1}, 
  PlotRange -> {Automatic, Automatic, {-0.1, 0.15}}, 
  PerformanceGoal -> "Quality"],
 {{t0, 0, "time"}, 0, 4, 0.01, Appearance -> "Labeled"},
 TrackedSymbols :> {t0}
 ]

NDSolve did give some warnings on console.

General::munfl: 6.17211*10^-153 7.58767*10^-160 is too small to represent 
as a normalized machine number; precision may be lost.

These look like from evaluating f at some points. May be with some additional options these can be eliminated.

For analytical solution, DSolve taking long time. So stopped it. It does not look like it can solve it analytically.

sol = DSolve[{pde, ic, bc}, u[x, y, t], {x, y, t}]

Mathematica 12.2 on windows 10

As an alternative, you can use HeatTransferPDEComponent and FEM to set up the model.

pde = HeatTransferPDEComponent[{u[x, y, t], t, {x, y}}, <||>] == 
   10*Sin[Pi*x]^10*Sin[Pi*y]^10*Sin[Pi*t]^10;
ic = u[x, y, 0] == 0;
bc = DirichletCondition[u[x, y, t] == 0, True];
sol = NDSolve[{pde, ic, bc}, 
   u, {x, y} \[Element] Rectangle[], {t, 0, 4}];

Also, note the Heat Transfer PDEs and Boundary Conditions guide page, the Heat Transfer tutorial, the Heat Transfer PDE models and the verification examples.