Spurious oscillations in solving diffusion problems using finite elements

Question

I've been struggling with this problem for a while so I hope someone can help me here.

I'm trying to solve the McNabb-Foster equations for hydrogen diffusion in metals using a simple 1D finite element simulation in Python. This works fine under isothermal conditions, but I can't make it work for anisothermic conditions (i.e. increasing temperature). The equations are essentially a non-linear Poisson equation: $$ \frac{dC}{dt}=D ∇^2 C-N\frac{dθ}{dt} $$ where $$ \frac{dθ}{dt}=κC(1-θ)-λθ $$

Where C is the hydrogen concentration, D, κ and λ are temperature-dependent terms, N is constant and θ is the so-called trap occupancy (traps are sites with high binding energy). Using the backward difference method the first equation can be expressed in matrix form as: $$ (\frac{1}{∆t}[M]+D[K])\{C\}_{t+∆t}={F}_{t+∆t}+\frac{1}{∆t}[M]\{C\}_t+[M]N\Bigl\{\frac{dθ}{dt}\Bigl\}_{t+∆t} $$ Where [M] and [K] are the system matrices. So far, so good, but the tricky part is how to incorporate θ. For the isothermal case I have done this again using the backward difference method, but taking θ from the previous iterative step: $$ \Bigl\{\frac{dθ}{dt}\Bigl\}_{t+∆t}=\frac{\{θ\}_{t+∆t}-\{θ\}_t}{Δt}=κ\{C\}_{t+∆t}\{1-θ\}_t-λ\{θ\}_t $$ Combining the above two equations gives: $$ (\frac{1}{∆t}[M]+D[K]+[M]Nκ\{1-θ\}_t )\{C\}_{t+∆t}=\{F\}_{t+∆t}+\frac{1}{∆t}[M]\{C\}_t+[M]Nλ\{θ\}_t $$ Below is the iterative part of my code where I have tried to implement this. Basically, I'm solving the equations twice; once to find the new concentrations C and a second time to find the flux F, which is actually what I need.

The trouble is that I get (spurious?) oscillations at the boundary for reasonable input parameters (see image). These disappear for very small time increments, but it makes the code too slow (I'm running it in a least-squares routine so it needs to be fast). The mesh size has no influence.

Does anyone have any suggestions how I can improve the stability?

for i in range(1, Iter+1):

     Clold = Cl # save Cl from last iteration
     T = Tstart + i * dt * dTdt # calculate temperature

     # calculate temperature dependent terms
     Dl = D0 * np.exp(-El / R / T)
     kappa = kappa0 * np.exp(-Et / R / T)
     lambd = lambda0 * np.exp(-Ed / R / T)

     # establish reverse difference equation
     LHS = sysm + dt * Dl * sysk + dt * Nt * kappa * sysm * (1-theta)
     RHS = np.dot(sysm, Clold) + dt * np.dot(sysm, (Nt * lambd * theta))
   
     # apply constraints & boundary conditions
     LHS[0, :] = 0.0
     LHS[0, 0] = 1.0
     RHS[0, 0] = 0.0

     # calculate new Cl
     Cl = np.linalg.solve(LHS, RHS)

     # calculate flux matrix   
     LHS = sysm + dt * Dl * sysk + dt * Nt * kappa * sysm * (1-theta)
     RHS = np.dot(sysm, Clold) + dt * np.dot(sysm, (Nt * lambd * theta))
     F = np.dot(LHS, Cl) - RHS
    
     # calculate new trap occupancies
     theta = theta + dt * (kappa * Cl * (1 - theta) - lambd * theta)

     # store flux for plotting
     Fstore = np.append(Fstore, np.array(-F[0] / dt))

EDIT: using an implicit scheme the last equation would look like this:

$$ (\frac{1}{∆t}[M]+D[K]+[M]Nκ\{1-θ\}_{t+∆t} )\{C\}_{t+∆t}=\{F\}_{t+∆t}+\frac{1}{∆t}[M]\{C\}_t+[M]Nλ\{θ\}_{t+∆t} $$

So now you have a second unknown term in the equation. The new θ can be given in terms of the new C, i.e.:

$$ θ_{t+∆t}=\frac{∆tκC_{t+∆t}+θ_t}{1+∆tκC_{t+∆t}+∆tλ} $$

but this makes the first equation non-linear. Any suggestion how to solve this?