I am wondering if there is an analytical solution for such a boundary value problem,
$$ \frac{\partial ^{2}c}{\partial x^2}+\frac{\partial c}{\partial x} -(u_{c1}(x)-u_{c2}(x))c =0 $$
where $u_{c1}(x)$ and $u_{c2}(x)$ are Heaviside function with the jump points at $x=c_1$ and $x=c_2$, respectively. The boundary condition is, $$ c(x=0)=c_0 $$ $$ c(x\rightarrow\infty)=0 $$ Any comments or suggestions are appreciated!
Let us follow the steps of this answer. For the sake of simplicity, let us assume $c_1<c_2$ and $c_0=0$. Using the definition of the Heaviside function, we can split the equation as follows: $$ \begin{cases} c''+c'=0, & 0\leq x<c_1\\ c''+c'-c=0, & c_1\leq x<c_2\\ c''+c'=0 & x\geq c_2. \end{cases}$$
Solving the first equation and applying the first boundary condition, $c(0)=0$, we get \begin{equation}\label{1} c(x)=\alpha_1(1-e^{-x}),\ 0\leq x<c_1.\tag{1} \end{equation}
In the same way, solving the $3^{rd}$ equation, and applying the second boundary condition, we have \begin{equation}\label{2} c(x)=\beta_3e^{-x}, \ x\geq c_2.\tag{2} \end{equation}
Now, we know that the second equation has the following general solution: \begin{equation}\label{3} c(x)=\alpha_2e^{\lambda_1x}+\beta_2e^{\lambda_2x}, \ c_1<x\leq c_2,\tag{3} \end{equation} where $\lambda_1=-\frac{1}{2} \left(1+\sqrt{5}\right)$ and $\lambda_2=\frac{1}{2} \left(-1+\sqrt{5}\right)$. If we seek for a continuous solution, we must evaluate the constant $c_1$ into \eqref{1} and \eqref{3} and then equals the result. The same for $c_2$. In carrying out these details, we obtain the following linear system: \begin{equation} \begin{cases} \alpha_2e^{\lambda_1c_1}+\beta_2e^{\lambda_2c_1}=\alpha_1(1-e^{-c_1})\\ \alpha_2e^{\lambda_1c_2}+\beta_2e^{\lambda_2c_2}=\beta_3e^{-c_2} \end{cases} \end{equation}
Solving this system, we will find $\alpha_2$ and $\beta_2$ as function of $\alpha_1$ and $\beta_3$, as follows: \begin{equation}\label{4}\alpha_2=-\frac{e^{-c_1-c_2} \left(-\beta_3 e^{c_1 +\lambda_2c_1}+\alpha_1 e^{c_1+c_2 \text{$\lambda_2$}+c_2}-\text{$\alpha_1 $} e^{c_2 \text{$\lambda_2 $}+c_2}\right)}{e^{c_1 \text{$\lambda_2 $}+c_2 \text{$\lambda_1 $}}-e^{c_1 \text{$\lambda_1 $}+c_2 \text{$\lambda_2 $}}}\tag{4}\end{equation} and \begin{equation}\label{5}\beta_2=\frac{\text{$\alpha_1 $} \left(1-e^{-c_1}\right) e^{c_2 \text{$\lambda_1 $}}-\text{$\beta_3 $} e^{c_1 \text{$\lambda_1 $}-c_2}}{e^{c_1 \text{$\lambda_2 $}+c_2 \text{$\lambda_1 $}}-e^{c_1 \text{$\lambda_1 $}+c_2 \text{$\lambda_2 $}}}.\tag{5}\end{equation}
So, the general solution to this problem will be given by $$c(x)=\begin{cases} \alpha_1(1-e^{-x}),& 0\leq x<c_1,\\ \alpha_2e^{\lambda_1x}+\beta_2e^{\lambda_2x},& \ c_1<x\leq c_2,\\ \beta_3e^{-x},& \ x\geq c_2, \end{cases} $$ where $\alpha_1$ and $\beta_3$ are arbitrary constants and $\alpha_2$ and $\beta_2$ are constants depending on $\alpha_2$ and $\beta_2$ given in \eqref{4} and \eqref{5}.
For instance, if we set $c_1=c_2=\alpha_1=\beta_3=1$, we will obtain the following graph: