Why is this Advection-convection model with insulating boundary losing mass?

Question

I am trying to model a 1-d advection-convection numerically, using an upwind scheme. I'm using the following equation to calculate the value of internal cells:

$$C_x^{t+1} = C_x^{t} + D\frac{\Delta t}{\Delta x^2} ( C_{x-1}^t-2 C_{x}^t+C_{x+1}^t) + u\frac{\Delta t}{\Delta x}(C_{x+1}^t-C_{x}^t)$$

and I'm using the Robin boundary condition (a.k.a. insulating boundary). In other words, I'm trying to guarantee no mass will leave the model environment. The boundary value is calculated by setting the leftmost and right most cells as:

$ C_1^{t+1} = D \frac{C_2^{t}}{u+D} $ and $ C_n^{t+1} = C_{n-1}^t\frac{(u\Delta x+D)}{D} $

Where $C$ is the concentration, $D$ the diffusion coefficient, $u$ is the velocity and $x$ is the current cell x-coordinate which goes from $1$ to $n$.

I understand this scheme should be numerically stable and accurate when both the following conditions are met:

$ u\frac{\Delta t}{\Delta x} + 2D\frac{\Delta t}{\Delta x^2} \leq1 $ and $\frac{u\Delta x}{D} < \frac{2}{1- u\frac{\Delta t}{\Delta x}} $

However, if I run this model with parameters $D=0.05$,$u=0.1$,$\Delta t=0.001$,$\Delta x=1$, there is some mass loss after 10,000 time steps, even tough the conditions above are met. It may seem a lot, but since my timesteps are very small it means less $t<10$ when that starts happening.

Any idea, of what could cause this loss of mass?

Below there's code to conduct this simulation in R.

transfun = function(time,state,parms) {
    D = parms$D
u = parms$u
    dt = parms$dt
dx = parms$dx
    n  = length(state)


    new = numeric(n)

    #Boundary cells
    state[1]= state[2]* D/(u+D)
    state[n]= state[n-1]* (u*dx+D)/D


    # Internal cells
    for(a  in 2:(n-1)) {
        state[a] = state[a] +D*dt/dx^2 * ( state[a-1] + -2*state[a] + state[a+1]) + 
                 u * dt/dx * (-state[a] + state[a+1])
    }

    #Boundary cells
    state[1]= state[2]* D/(u+D)
    state[n]= state[n-1]* (u*dx+D)/D


    return(list(state))
}

## Hand made loop 
# Ensure parameters follow the stability rule: 0<= C^2 <= 2s <= 1 
# where C = u*dt/dx and s = D*dt/dx^2
D=0.05
u=0.1
dt=0.001
dx=1
C = u*dt/dx
s = D*dt/(dx^2)
n = 15
T = 10000
print(C+2*s < 1)
print(u*dx/D<2/(1-C))



resul = matrix(NA,T,n)
resul[1,] = rep(0,n)
resul[1,ceiling(15/2)] = 1
for(a in seq(2,T,1)) {
    resul[a,] = transfun(1,state=resul[a-1,],parms=list(D=D,u=u,dt=dt,dx=dx))[[1]]
}

,