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.
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!