先上问题

res3 = NDSolveValue[{Laplacian[u[x, y, t], {t}] ==
Laplacian[u[x, y, t], {x, y}],
DirichletCondition[u[x, y, t] == p[x, y], t == 0],
D[u[x, y, t], t] == NeumannValue[0, t == 0],
DirichletCondition[D[u[x, y, t], t] == 0, t == 0],
DirichletCondition[u[x, y, t] == p[x, y], x == 3],
DirichletCondition[u[x, y, t] == p[x, y], x == -3],
DirichletCondition[u[x, y, t] == p[x, y], y == 3],
DirichletCondition[u[x, y, t] == p[x, y], y == -3]},
u, {t, 0, 1}, {x, -3, 3}, {y, -3, 3}]

其实吧，是你的代码写错了。。。另外下一次记得把p[x,y]的代码也贴上。。。
p[x_, y_] :=
3 (1 - x)^2 Exp[-x^2 - (y + 1)^2] -
10 (x/5 - x^3 - y^5) Exp[-x^2 - y^2] - 1/3 Exp[-(x + 1)^2 - y^2]
res3 = NDSolve[{Laplacian[u[x, y, t], {t}] -
Laplacian[u[x, y, t], {x, y}] == 0,
DirichletCondition[u[x, y, t] == p[x, y], True],
u[x, y, 0] == p[x, y], Derivative[0, 0, 1][u][x, y, 0] == 0},
u, {t, 0, 1}, {x, -3, 3}, {y, -3, 3}, Method -> {
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.1},
"IntegrationOrder" -> 5}}}]
实际上Method那个代码完全可以删去得到结果，加上是指定使用FEM方法得到更好的效果，具体见文档
FEMDocumentation/tutorial/FiniteElementProgramming#752095140