On the analytical solution of a linear system of coupled partial differential equations with Dirac functions using Fourier series expansion method I am trying to find an analytical solution of the following system of coupled partial differential equations:
\begin{align}
  \rho_{,t} &= 
   \rho_{,xx} - \frac{\phi_{,x}}{2} + P(x) \, , \\
  \phi_{,t} &=
 \rho_{,x} - \phi\, ,
\end{align}
where comma in indices stads for a partial derivative. In addition, 
$$
P(x) =  \delta(x) 
 - \frac{1}{2} \Big( \delta\left(x-\tfrac{1}{2}\right)
 +\delta\left(x+\tfrac{1}{2}\right)
  \Big)  \, .
$$
This system is subject to the boundary conditions $\rho(x= \pm \frac{1}{2}, t)=0$ (zero displacement at boundaries), in addition to the initial conditions $\rho(x,t=0) = \phi(x,t=0)=0$.

For that purpose, I tried to expand the fields in Fourier series in such a way the BC's above are satisfied. 
  Accordingly,
  \begin{align}
  \rho(x,t) &= \sum_{n \ge 1} \rho_n(t) \left( 1 + c_n(x) \right) \, , \\
  \phi(x,t) &= \sum_{n \ge 1} \phi_n(t) \, s_n(x) \, , 
 \end{align}
  where $\rho_n$ and $\phi_n$ are unknown time-dependent series coefficients.
  Moreover,
      \begin{align}
  c_n(x) &= \cos \left( 2(2n-1)\pi x \right) \, , \\
  s_n(x) &= \sin \left( 2(2n-1)\pi x \right) \, .
 \end{align}

By making use of the Fourier series expansion of the delta Dirac function,
\begin{equation}
 \delta(x)=1
 + 2 \sum_{n\ge 1} \cos \left( 2\pi n x \right) \, .
\end{equation}
it is easy to show that
$$
P(x) = 4\sum_{n\ge 1} c_n(x) \, .
$$
Upon substitution in the original equations, we obtain
\begin{align}
    \sum_{n\ge 1} \frac{\mathrm{d} \rho_n}{\mathrm{d} t} \big( 1 + c_n(x) \big)
   &= 
   \sum_{n\ge 1} \Bigg( -H_n \left( H_n \, \rho_n 
   +\frac{\phi_n}{2} \right) + 4
    \Bigg) c_n(x) \, , \tag{1} \\
 \sum_{n\ge 1} \frac{\mathrm{d} \phi_n}{\mathrm{d} t} \, s_n(x) 
 &= - \sum_{n\ge 1} \left( H_n\, \rho_n + \phi_n \right) s_n (x) \, , \tag{2}
\end{align}
where we have defined $H_n = 2(2n-1)\pi$.
By multiplying both members of Eqs. (1) and (2) by $c_m(x)$ and $s_m(x)$, respectively, integrating with respect to $x$ over $-1/2$ and $1/2$, and making use of the orthogonality relations, 
\begin{equation}
 \int_{-\frac{1}{2}}^{\frac{1}{2}} c_n(x) c_m(x) \, \mathrm{d} x 
 =\int_{-\frac{1}{2}}^{\frac{1}{2}} s_n(x) s_m(x) \, \mathrm{d} x 
 = \frac{1}{2} \, \delta_{mn} , 
\end{equation}
only the terms $m=n$ of the infinite sums remain, to obtain
    \begin{align}
   \frac{\mathrm{d} \rho_n}{\mathrm{d} t} 
  &= - H_n \bigg( H_n \, \rho_n
  + \frac{\phi_n}{2} \bigg) + 4   \, , \tag{3} \\
   \frac{\mathrm{d} \phi_n}{\mathrm{d} t} 
  &= - H_n \, \rho_n - \phi_n \, . \tag{4}
 \end{align}
The system composed of Eqs. (3) and (4) can readily be solved using Laplace transforms, or alternatively, using Maple or Mathematica, to obtain
\begin{align}
 \rho_n(t) &= \frac{8}{H_n^2}
 \Bigg( 1+\frac{\delta_n^+}{T_n} \left( \delta_n^--1 \right)e^{-\delta_n^-t}
 -\frac{\delta_n^-}{T_n} \left( \delta_n^+ -1\right)e^{-\delta_n^+ t} \Bigg) \, , \\
 \phi_n (t) &= -\frac{8}{H_n} 
 \left( 1 -\frac{\delta_n^+}{T_n} \, e^{-\delta_n^- t} + \frac{\delta_n^-}{T_n} \, e^{-\delta_n^+ t} \right) \, , 
\end{align}
where we have defined
    \begin{align}
  T_n &= \sqrt{1+H_n^4} \, , \\
  \delta_n^\pm &= \frac{1}{2} \left( H_n^2+1 \pm T_n \right) \, .
 \end{align}


The problem one is facing is that this solution does not satisfy Eqs. (1) when evaluating the left and right-hand sides at some arbitrary points, say $x=1/4$ and $t=1/10$.
    Eq. (2) on the other hand is well satisfied. 


I would be grateful if someone here could find some time to have a look at my calculations and let me know whether some math steps/ assumptions, notably the ones leading to Eqs. (3) and (4), are wrong. 
Any ideas/ suggestions are most welcome -- Thank you

Clearly, $\delta_n^\pm \ge 0$, and thus both $\rho_n$ and $\phi_n$ are well behaved in the limit as $t\to\infty$, to obtain
    \begin{align}
  \lim_{t\to\infty} \rho(x,t) &=
  \frac{2}{\pi^2} \sum_{n\ge 1} \frac{1+c_n(x)}{(2n-1)^2} \, , \\
  \lim_{t\to\infty} \phi(x,t) &= -\frac{4}{\pi} \sum_{n\ge 1} 
  \frac{s_n(x)}{2n-1} \, , 
 \end{align}
which correspond, respectively, to the Fourier series representation of the triangle and square waves functions of frequency $2\pi$.
The method seems to work quite nicely in the steady limit where the solution can easily be obtained and is given by
    \begin{align}
  \rho(x) &= P_0 \left( \frac{1}{2} - |x| \right) \, , \\
  \phi(x) &= -P_0 \, \mathrm{sgn} (x) \, ,
 \end{align}
where $\mathrm{sgn}(x) := x/|x|$ denotes the sign function.
 A: the solution of your problem can more conveniently by obtained using Fourier transoms in space, followed by Laplace transform of the temporal variable to solve the resulting ODE. 
    \begin{align}
  \hat{\rho}_{,t} &= -q^2 \hat{\rho} - \frac{iq}{2} \, \hat{\phi} 
  + 2 \sin \left( \frac{q}{4} \right)^2 
   \, ,  \\
  \hat{\phi}_{,t} &= iq \hat{\rho} - \hat{\phi} \,.
 \end{align}
Employing the initial conditions $\hat{\rho}(q,t=0)=\hat{\phi}(q,t=0)=0$, we obtain
    \begin{align}
  s \, \hat{\rho}(q,s) &= -q^2 \hat{\rho}(q,s) - \frac{iq}{2} \, \hat{\phi} (q,s)
    + \frac{2}{s} \,  \sin \left( \frac{q}{4} \right)^2  \, , \\
  s \, \hat{\phi}(q,s) &= iq \hat{\rho}(q,s) - \hat{\phi}(q,s) \, . 
 \end{align}
Solving for $\hat{\rho} (s)$ and $\hat{\phi} (s)$ yields
\begin{equation}
 \hat{\rho} (q, s) = \frac{4(1+s)}{Q} \, \sin \left( \frac{q}{4} \right)^2  \, , \qquad
 \hat{\phi} (q, s) = \frac{4iq}{Q} \, \sin \left( \frac{q}{4} \right)^2  \, , 
\end{equation}
where we have defined
\begin{equation}
 Q = s \big( 2s^2 + 2(1+ q^2)s + q^2 \big) \, .
\end{equation}
By applying the inverse Laplace transform, we obtain
    \begin{align}
  \hat{\rho} (q,t) &= \frac{4}{q^2} \, \sin \left( \frac{q}{4} \right)^2 \bigg[  
  1-e^{-\beta t}
  \left( \cosh \left( \alpha t \right) 
  + \frac{1}{2\alpha} \, \sinh \left( \alpha t\right) 
   \right)
  \bigg] \, , \notag \\
  \hat{\phi} (q,t) &= \frac{4i }{q} \, \sin \left( \frac{q}{4} \right)^2 \left[
  1- e^{-\beta t } 
  \left( \cosh(\alpha t) + \frac{\beta}{\alpha} \, \sinh (\alpha t) \right)
  \right] \, , 
 \end{align}
where we have defined
\begin{equation}
 \alpha = \frac{1}{2} \sqrt{1+q^4} \, , \qquad
 \beta = \frac{1}{2} \left( 1+q^2 \right) \, .
\end{equation}
The solution in real space can readily be obtained (via numerical integration) upon inverse Fourier transformation.
Hope this could help a bit.
Why your method solution, even though it look correct, does not work still is an open question!
