On the "Derivation of the Electromagnetic Lagrangian density" In the most upvoted answer here : Deriving Lagrangian density for electromagnetic field, how do we know that  equations (015) and (016) therein
\begin{equation}
\boxed{\:
\dfrac{\partial }{\partial t}\left(\nabla \boldsymbol{\cdot} \mathbf{A}\right)+\nabla \boldsymbol{\cdot} \left(\boldsymbol{\nabla}\phi \right)  -\left(-\frac{\rho}{\epsilon_{0}}\right) =0  
\:}
\tag{015} 
\end{equation}
\begin{equation}
\boxed{\:
\dfrac{\partial}{\partial t}\left(\frac{\partial \mathrm{A}_{k}}{\partial t}+\frac{\partial \phi}{\partial x_{k}}\right)+\nabla \boldsymbol{\cdot}  \left[c^{2}\left(\frac{\partial \mathbf{A}}{\partial x_{k}}- \boldsymbol{\nabla}\mathrm{A}_{k}\right)\right] -\frac{\mathrm{j}_{k}}{\epsilon_{0}}=0
 \:}
\tag{016} 
\end{equation}
will give the Lagrangian density and not the Lagrangian function?
 A: The Lagrange function is used for mechanical systems with finite or at most denumerably infinite number of degrees of freedom
\begin{equation}
L\left(q_{\boldsymbol i},\overset{\,\centerdot}{q}_{\boldsymbol i},t\right) \qquad i=1,2,3,\cdots
\tag{01}\label{01}     
\end{equation}
Each value of the discrete index $\;i\;$ corresponds to a different one of the generalized coordinates $\;q_{\boldsymbol i}$.
The Euler-Lagrange equations (equations of motion) for such systems are expressed via the function \eqref{01}
\begin{equation}
\dfrac{\mathrm d}{\mathrm dt}\left(\dfrac{\partial L}{\partial\overset{\,\centerdot}{q}_{\boldsymbol i} }\right)-\dfrac{\partial L}{\partial q_{\boldsymbol i}}\boldsymbol{=}0 \qquad i=1,2,3,\cdots
\tag{02}\label{02}     
\end{equation}
There are some mechanical problems, however, that involve continuous systems, as, for example, the problem of a vibrating continuous elastic rod. Here each point of the continuous rod partakes in the longitudinal oscillations, and the complete motion can only be described by specifying the position coordinates of all points. In this continuous one dimensional case the role of the discrete generalized coordinates $\;q_{\boldsymbol i}\;$ is played by a real scalar function $\;\eta\left(x\right)$, the displacement along the rod axis of a point with position coordinate $\;x$. This position coordinate $\;x\;$ is an index of the generalized coordinates $\;\eta\left(x\right)\;$ in a way similar to the discrete case where $\;i\;$ is an index of the generalized coordinates $\;q_{\boldsymbol i}$. Since $\;\eta\;$ depends also upon the continuous variable $\;t$, we should perhaps write more accurately $\;\eta\left(x,t\right)$, indicating that $\;x$, like $\;t$, can be considered as a parameter entering into the Lagrangian.        
The Lagrangian density for any one-dimensional continuous system, like the rod above, would appear as a function of the form
\begin{equation}
\mathcal{L}\left(\eta,\dfrac{\mathrm d \eta}{\mathrm dx},\dfrac{\mathrm d \eta}{\mathrm dt},x,t\right)
\tag{03}\label{03}    
\end{equation}
while the Lagrange function  $\;L\;$ as the integral of $\;\mathcal{L}\;$ over $\;x$
\begin{equation}
L\boldsymbol{=}\int \mathcal{L}\mathrm dx
\tag{04}\label{04}    
\end{equation}
Now, Hamilton's principle for a continuous system is the same 
\begin{equation}
\delta I=\delta \int\limits_{\bf 1}^{\bf 2} L\mathrm dt\boldsymbol{=}\delta  \int\limits_{\bf 1}^{\bf 2}\!\!\int\mathcal{L}\mathrm dx \mathrm dt\boldsymbol{=}0
\tag{05}\label{05}    
\end{equation}
After elaborating  the Hamilton's principle could be written as
\begin{equation}
\int\limits_{\boldsymbol t_{1}}^{\boldsymbol t_{2}}\int\limits_{\boldsymbol x_{1}}^{\boldsymbol x_{2}}\mathrm dx \mathrm dt \left[\dfrac{\mathrm d}{\mathrm dt}\left(\dfrac{\partial\mathcal{L}  }{\partial\tfrac{\mathrm d\eta}{\mathrm dt} }\right)\boldsymbol{+}\dfrac{\mathrm d}{\mathrm dx}\left(\dfrac{\partial\mathcal{L}}{\partial\tfrac{\mathrm d\eta}{\mathrm dx}}\right)\boldsymbol{-}\dfrac{\partial\mathcal{L}}{\partial\eta}\right]\left(\dfrac{\partial\eta}{\partial \alpha}\right)_{0}\boldsymbol{=}0
\tag{06}\label{06}    
\end{equation}
where $\;\alpha\;$ the parameter of a family of $\;\eta\;$ functions of the form
\begin{equation}
\eta\left(x,t;\alpha\right)\boldsymbol{=}\eta\left(x,t;0\right)\boldsymbol{+}\alpha\zeta\left(x,t;0\right)
\tag{07}\label{07}    
\end{equation}
Here $\;\eta\left(x,t;0\right)\;$ stands for the correct function that will satisfy Hamilton's principle,
and $\;\zeta\left(x,t\right)\;$ is any well-behaved function that vanishes at the end points in $\;t\;$ and  $\;x$. The arbitrary nature of the varied path implies the vanishing of the expression in the brackets under the integral in \eqref{06}
\begin{equation}
\dfrac{\mathrm d}{\mathrm dt}\left(\dfrac{\partial\mathcal{L}  }{\partial\tfrac{\mathrm d\eta}{\mathrm dt} }\right)\boldsymbol{+}\dfrac{\mathrm d}{\mathrm dx}\left(\dfrac{\partial\mathcal{L}}{\partial\tfrac{\mathrm d\eta}{\mathrm dx}}\right)\boldsymbol{-}\dfrac{\partial\mathcal{L}}{\partial\eta}\boldsymbol{=}0
\tag{08}\label{08}    
\end{equation}
So

In continuous systems the Euler-Lagrange equations or equations of motion are expressed by the Lagrangian Density  $\;\mathcal{L}\;$ and not by the Lagrange function $\;L$.

In a 3-dimensional case the Lagrangian Density  $\;\mathcal{L}\;$ has the form
\begin{equation}
\mathcal{L}\left(\eta,\overset{\,\centerdot}{\eta},\boldsymbol{\nabla}\eta,\mathbf{x},t\right)
\tag{09}\label{09}    
\end{equation}
and the Euler-Lagrange equations are
\begin{equation}
\dfrac{\partial }{\partial t}\left(\dfrac{\partial \mathcal{L}}{\partial \overset{\,\centerdot}{\eta}}\right)\boldsymbol{+}\boldsymbol{\nabla\cdot}  \left[\dfrac{\partial \mathcal{L}}{\partial\left(\boldsymbol{\nabla}\eta\right)}\right]\boldsymbol{-}\frac{\partial \mathcal{L}}{\partial\eta}\boldsymbol{=}0
\tag{10}\label{10}    
\end{equation}


Reference : $^{\prime\prime}$Classical Mechanics$^{\prime\prime}$ by H.Goldstein-C.Poole-J.Safko, 3rd Edition 2000, $\S$Chapter 13 : Introduction to the Lagrangian and Hamiltonian Formulations for Continuous Systems and Fields.

