I have this equation for $\rho$ $$ \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial x}(\rho u) = \varepsilon\frac{\partial^2 \rho}{\partial x^2},\quad 0\leq x\leq 1,\quad 0\leq t\leq T\\ \frac{\partial \rho}{\partial x}\Bigg|_{x=0;1} = 0\\ \rho|_{t=0} = \varphi(x)\\ u = u(x,t) $$ I'm trying to get a numerical answer with implicit finite difference methods. I tried to build a solution for this, but I can't get it to work. Here is what I tried. $$ x_0<\ldots<x_M,\quad x_i=ih,\quad x_M=Mx=1\\ t_0<\ldots<t_N,\quad t_j=j\tau,\quad t_N=Nt=T\\ \frac{\rho_{i}^{j+1}-\rho_{i}^{j}}{\tau} + u_i^{j}\frac{\rho_{i+1}^{j+1} - \rho_{i}^{j+1}}{h} + \rho_{i}^{j}\frac{u_{i+1}^{j+1} - u_{i}^{j+1}}{h} = \varepsilon\frac{\rho_{i+1}^{j+1} - 2\rho_{i}^{j+1} + \rho_{i-1}^{j+1}}{h^2}\\ \frac{\rho_1^{j+1} - \rho_0^{j+1}}{h}=0,\quad \frac{\rho_M^{j+1} - \rho_{M-1}^{j+1}}{h}=0 $$ And then I tried to solve it programaticaly using tridiagonal matrix algorithm. First, I compose it to $$ A_i^{j+1}\rho_{i-1}^{j+1} - C_i^{j+1}\rho_i^{j+1} + B_i^{j+1}\rho_{i+1}^{j+1} = -F_i^{j+1} $$ and find coefficients $$ A_i^{j+1} = -\frac{\varepsilon}{h^2},\\ C_i^{j+1} = -\frac{1}{\tau} + \frac{u_i^j}{h} - \frac{2\varepsilon}{h^2},\\ B_i^{j+1} = \frac{u_i^j}{h} - \frac{\varepsilon}{h^2},\\ F_i^{j+1} = -\frac{\rho_i^j}{\tau} + \rho_i^j\frac{u_{i+1}^{j+1} - u_i^{j+1}}{h}. $$ Let $$ \rho_i^{j+1}=\alpha_{i}^{j+1}\rho_{i+1}^{j+1} + \beta_i^{j+1} $$ and by substitution I get $$ \alpha_i^{j+1} = \frac{B_i^{j+1}}{C_i^{j+1} - \alpha_{i-1}^{j+1}A_{i}^{j+1}},\\ \beta_i^{j+1} = \frac{A_i^{j+1}\beta_{i-1}^{j+1} + F_i^{j+1}}{C_i^{j+1} - \alpha_{i-1}^{j+1}A_{i}^{j+1}} $$ And then from the first boundary condition I found that $\alpha_0^{j} = 1$ and $\beta_0^{j} = 0$. And from second boundary condition I got $\rho_M^{j} = \frac{\beta_{M-1}^{j}}{1 - \alpha_{M-1}^{j}}$.
But it doesn't work and gives values that are impossible for given situation (I tested by adding f(x,t) to the equation and making $\rho = e^{-t}\cos(\pi x)$. And the numerical solution behaves completely different. I think, that there is a problem with how I approximated the boundary conditions (specificaly how I applied the second one) or the partial differential of $\rho u$. Can someone suggest me a stable solution for this problem?
Here is what I get with this solution: solution. First plot is what I get, the second one - actual solution.
I've read so much literature on this topic for the past couple of weeks, but found nothing, that could help me. Also, I'm sorry if my wording isn't good, I'm not native to english.