Sorry, my patch above is not quite right, you of course need to use the copied A for the solution too:
--- a/main.jl
+++ b/main.jl
@@ -135,8 +135,9 @@ for t in 1:Δt:T
update!(ch, t) # load current dbc values from input_exp
global b = Δt*f + M * c_n # get the discrete rhs
- apply!(A, b, ch) # apply time-dependent dbc
- global c = A \ b # solve the current time step
+ copyA = copy(A)
+ apply!(copyA, b, ch) # apply time-dependent dbc
+ global c = copyA \ b # solve the current time step
push!(store, c) # store current solution
global c_n = copy(c) # update time step
which results in this:
