We investigate the inverse problem of numerically identifying unknown initial temperatures in a heat equation with dynamic boundary conditions whenever some overdetermination data is provided after a final time. This is a backward parabolic problem which is severely ill-posed. As a first step, the problem is reformulated as an optimization problem with an associated cost functional. Using the weak solution approach, an explicit formula for the Fr\'echet gradient of the cost functional is derived from the corresponding sensitivity and adjoint problems. Then the Lipschitz continuity of the gradient is proved. Next, further spectral properties of the input-output operator are established. Finally, the numerical results for noisy measured data are performed using the regularization framework and the conjugate gradient method. We consider both one- and two-dimensional numerical experiments using finite difference discretization to illustrate the efficiency of the designed algorithm. Aside from dealing with a time derivative on the boundary, the presence of a boundary diffusion makes the analysis more complicated. This issue is handled in the 2-D case by considering the polar coordinate system. The presented method implies fast numerical results.