$$
\newcommand{\uex}{{u_{\small\mbox{e}}}}
\newcommand{\half}{\frac{1}{2}}
\newcommand{\halfi}{{1/2}}
\newcommand{\xpoint}{\boldsymbol{x}}
\newcommand{\normalvec}{\boldsymbol{n}}
\newcommand{\Oof}[1]{\mathcal{O}(#1)}
\newcommand{\Ix}{\mathcal{I}_x}
\newcommand{\Iy}{\mathcal{I}_y}
\newcommand{\It}{\mathcal{I}_t}
\newcommand{\setb}[1]{#1^0} % set begin
\newcommand{\sete}[1]{#1^{-1}} % set end
\newcommand{\setl}[1]{#1^-}
\newcommand{\setr}[1]{#1^+}
\newcommand{\seti}[1]{#1^i}
\newcommand{\Real}{\mathbb{R}}
$$
Discretization
$$
[D_tD_t u = c^2(D_xD_x u + D_yD_yu) + f]^n_{i,j,k},
$$
Written out in detail:
$$
\begin{align*}
\frac{u^{n+1}_{i,j} - 2u^{n}_{i,j} + u^{n-1}_{i,j}}{\Delta t^2}
&= c^2
\frac{u^{n}_{i+1,j} - 2u^{n}_{i,j} + u^{n}_{i-1,j}}{\Delta x^2}
+ \nonumber\\
&\quad c^2\frac{u^{n}_{i,j+1} - 2u^{n}_{i,j} + u^{n}_{i,j-1}}{\Delta y^2}
+ f^n_{i,j},
\end{align*}
$$
\( u^{n-1}_{i,j} \) and \( u^n_{i,j} \) are known, solve for
\( u^{n+1}_{i,j} \):
$$ u^{n+1}_{i,j} = 2u^n_{i,j} + u^{n-1}_{i,j} + c^2\Delta t^2[D_xD_x u + D_yD_y u]^n_{i,j}$$
