Dtyjlui

I'm trying to obtain an approximation of the solution of the following equation:
$$
leftlbrace 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 + jDelta x,hspace{.2cm} j in lbrace0,1,...,N+1rbrace \
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}}{2Delta 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)}{2Delta x} right) + u_j(t)left( beta_x(t,xj) - 2frac{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)}{2Delta x} right)
$$
Naming $a_1(t,x_j) = frac{alpha(t,x_j)}{(Delta x)^2} + frac{beta(t,x_j)}{2Delta x}$, $a_2(t,x_j) = beta_x(t,xj) - 2frac{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)}{2Delta 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 $Ntimes 1$. A is the following $Ntimes 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 + iDelta t, hspace{.2cm} i in lbrace 0, 1, ..., Mrbrace \
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)}{2Delta x} right) + u_j(t)left( beta_x(t,xj) - 2frac{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)}{2Delta 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)}{2Delta x} right) + \ u_j(t)left( beta_x(t,xj) - 2frac{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)}{2Delta 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.