I'm trying to obtain an approximation of the solution of the following equation: $$ \left\lbrace \begin{array}{l,l} u_t = \alpha u_{xx} + (\beta u)_x, & u,\alpha ,\beta \in [T_0,T_f]\times [X_0,X_f]\\ u(t=T_0,x) = u_0(x), & \forall x \in [X_0,X_f] \\ u(t,x=X_0) = g_0(t), u(t,x=X_f) = g_f(t), & \forall t \in [T_0,Tf] \end{array}\right. $$ With $u_t = \frac{\partial u}{\partial t}$, $(\beta u)_x = \frac{\partial (\beta u)}{\partial x}$ and $u_{xx} = \frac{\partial^2 u}{(\partial x)^2}$.
The thing is, there must be something wrong since whenever $\beta$ is negative the approximations end up blowing up several times (so far it goes well for as long as this is avoided). I'll explain what I used for obtaining this approximations in case you can see if there's something I missed or simply is mistaken.
I used a finite difference method, by first applying the method of lines leaving the $t$ variable continuous as $x$ is discretized with $N+2$ nodes as follows: $$ \begin{array}{l} \Delta x = \frac{X_f - X_0}{N+1} \\ x_j = X_0 + j\Delta x,\hspace{.2cm} j \in \lbrace0,1,...,N+1\rbrace \\ u_j(t) \simeq u(t,x_j) \end{array} $$ Using central difference for $u_{xx}(x)$ ($u_{xx}\simeq \frac{u(x+\Delta x) - 2u(x) + u(x-\Delta x)}{(\Delta x)^2}$) the equation results as follows: $$ \frac{\partial u_j}{\partial t}(t) = \alpha(t,x_j)\frac{u_{j+1} - 2u_j + u_{j-1}}{(\Delta x)^2}(t) + \beta(t,x_j)\frac{u_{j+1}-u_{j-1}}{2\Delta x}(t) + \beta_x(t,x_j) u_j(t) = \\ = u_{j-1}(t)\left( \frac{\alpha(t,x_j)}{(\Delta x)^2} - \frac{\beta(t,x_j)}{2\Delta x} \right) + u_j(t)\left( \beta_x(t,xj) - 2\frac{\alpha(t,x_j)}{(\Delta x)^2} \right) + u_{j+1}(t)\left( \frac{\alpha(t,x_j)}{(\Delta x)^2} +\frac{\beta(t,x_j)}{2\Delta x} \right) $$ Naming $a_1(t,x_j) = \frac{\alpha(t,x_j)}{(\Delta x)^2} + \frac{\beta(t,x_j)}{2\Delta x}$, $a_2(t,x_j) = \beta_x(t,xj) - 2\frac{\alpha(t,x_j)}{(\Delta x)^2}$ and $a_3(t,x_j) = \frac{\alpha(t,x_j)}{(\Delta x)^2} +\frac{\beta(t,x_j)}{2\Delta x}$ the following system of equations is obtained: $$ U_t(t) = A(t)U(t) $$ With: $$ \begin{array}{l} U_t(t) \simeq (u_t(t,x_1),...,u_t(t,x_N))' \\ U = (u_1(t),...,u_N(t))' \\ g = (a_1(t,x_1)g_0(t),0,0,...,0,a_3(t,x_N)g_f(t))' \end{array} $$ All of them being $N\times 1$. A is the following $N\times N$ matrix:
$$ \left(\begin{array}{c,c,c,c,c,c,c,c} a_2(t,x1) & a_3(t,x_1) & 0 & 0 & ... & ...& ... & 0 \\ a_1(t,x_2) & a_2(t,x_2) & a_3(t,x_2) & 0 & ... & ... & ... & 0 \\ 0 & a_1(t,x_3) & a_2(t,x_3) & a_3(t,x_3) & 0 & ... & ... & 0 \\ 0 & 0 & \ddots & \ddots & \ddots & 0 & ... & 0 \\ \vdots & 0 & 0 & \ddots & \ddots & \ddots & ... & \vdots \\ 0 & ... & ... & 0 & 0 & a_1(t,x_{N-1}) & a_2(t,x_{N-1}) & a_3(t,x_{N-1}) \\ 0 & 0 & ... & ... & ... & 0 & a_1(t,x_N) & a_2(t,x_N) \end{array}\right) $$ Then I discretized time with $M+1$ nodes as follows: $$ \begin{array}{l} \Delta t = \frac{T_f - T_0}{M}\\ t_i = T_0 + i\Delta t, \hspace{.2cm} i \in \lbrace 0, 1, ..., M\rbrace \\ u_j^i \simeq u(t_i,x_j) \end{array} $$ Applying forward Euler ($\frac{\partial u_j^i}{\partial t} \simeq \frac{u_j^{i+1} - u_j^i}{\Delta t}$) to the equation obtained after space discretization: $$ \frac{u_j^{i+1} - u_j^i}{\Delta t} = u_{j-1}(t)\left( \frac{\alpha(t,x_j)}{(\Delta x)^2} - \frac{\beta(t,x_j)}{2\Delta x} \right) + u_j(t)\left( \beta_x(t,xj) - 2\frac{\alpha(t,x_j)}{(\Delta x)^2} \right) + u_{j+1}(t)\left( \frac{\alpha(t,x_j)}{(\Delta x)^2} +\frac{\beta(t,x_j)}{2\Delta x} \right) $$ Hence, we can state $u_j^{i+1}$ as: $$ u_j^{i+1} = u_j^i + \Delta t \left(u_{j-1}(t)\left( \frac{\alpha(t,x_j)}{(\Delta x)^2} - \frac{\beta(t,x_j)}{2\Delta x} \right) + \\ u_j(t)\left( \beta_x(t,xj) - 2\frac{\alpha(t,x_j)}{(\Delta x)^2} \right) + u_{j+1}(t)\left( \frac{\alpha(t,x_j)}{(\Delta x)^2} +\frac{\beta(t,x_j)}{2\Delta x} \right)\right) $$ With the following notation: $$ \begin{array}{l} U^i = (u_1^i,u_2^i,...,u_N^i)' \\ h^i = (a_1(t_i,x_ 1)g_0(t_i)\Delta t,0,...,0,a_3(t_i,x_N)g_f(t_i)\Delta t )' \end{array} $$ $U^i$ and $h^i$ with dimensions $N \times 1$. Considering $I_N$ the $N \times N$ identity matrix we have: $$ U^{i+1} = \left(I_N + \Delta t A(t_i)\right)U^i + h^i, \hspace{.2cm} 1 \leq i \leq N $$
I plotted the approximations at some time values and, as I said before, whenever $\beta$ has negative values several blow ups appear. I would like to know if there's something wrong with the method I used or mistakes during the process so I can determine whether the error lies in the method or the code.
Thank you in advance!
P.S: The discretizations are taken in a way the Courant-Friedrichs-Levy condition is satisfied, that is not the issue.