diff --git a/examples/bruss1d.jl b/examples/bruss1d.jl
new file mode 100644
index 000000000..4df7035eb
--- /dev/null
+++ b/examples/bruss1d.jl
@@ -0,0 +1,202 @@
+using ODE,DASSL
+
+const T = Float64 # the datatype used, probably Float64
+Tarr = SparseMatrixCSC
+const N = 500 # active gridpoints
+const dof = 2N # degrees of freedom (boundary points are not free)
+dx = 1/(N+1)
+x = dx:dx:1-dx
+dae = 0  # index of DAE, ==0 for ODE
+
+# parameters
+alpha = 1/50
+const gamma = alpha/dx^2
+
+# BC
+const one_ = one(T)
+const three = 3*one_
+ubc1() = one_
+ubc2() = one_
+vbc1() = three
+vbc2() = three
+#IC
+u0(x) = 1 + 0.5*sin(2π*x)
+v0(x) = 3 * ones(length(x))
+
+inds(i) = (2i-1, 2i)
+function fn!(res, y, dydt)
+    # the ode function res = f(y,dydt)
+    #
+    # Note the u and v are staggered:
+    # u = y[1:2:end], v=y[2:2:end]
+    
+    # setup for first step
+    i = 1
+    iu, iv = inds(i)
+    ui0 = ubc1(); vi0 = vbc1() # BC
+    ui1 = y[iu];   vi1 = y[iv]
+    ui2 = y[iu+2]; vi2 = y[iv+2]
+    while i<=N
+        res[iu] = 1 + ui1^2*vi1 - 4*ui1 + gamma*(ui0 - 2*ui1 + ui2) - dydt[iu]
+        res[iv] = 3*ui1 - ui1^2*vi1     + gamma*(vi0 - 2*vi1 + vi2) - dydt[iv]
+        # setup for next step:
+        i += 1
+        iu, iv = inds(i)
+        ui0 = ui1; vi0 = vi1
+        ui1 = ui2; vi1 = vi2
+        if i<N
+            ui2 = y[iu+2]
+            vi2 = y[iv+2]
+        else # on last step do BC
+            ui2 = ubc2()
+            vi2 = vbc2()
+        end
+    end
+    return nothing
+end
+fn!(;T_::Type=T, dof_=dof) = zeros(T_,dof_)
+function fn(t, y, ydot)
+    res = fn!()
+    fn!(res, y, ydot)
+    return res
+end
+function fn_ex(t, y)
+    res = fn!()
+    fn!(res, y, 0*y)
+    return res
+end
+
+
+
+# Jacobian is a banded matrix with upper and lower bandwidth 2,
+# c.f. W&H p.148, but Julia does not support banded matrices yet.
+# Thus use a sparse matrix.
+
+function jac!(dfdy, y, ydot, a)
+    # The Jacobian of fn: = df/dy + a* df/dydot
+    ii = 1 # direct index into dfdy.nzval
+    for i=1:N # iterate over double-columns
+        iu, iv = inds(i)
+        ui = y[iu]
+        vi = y[iv]
+        # iu: column belonging to d/du(i)
+        if iu-2>=1
+            # du'(i-1)/du(i)
+            dfdy.nzval[ii] = gamma
+            ii +=1
+        end
+        if iv-2>=1
+            # dv'(i-1)/du(i)==0
+            dfdy.nzval[ii] = 0
+            ii +=1
+        end
+        # du'(i)/du(i)
+        dfdy.nzval[ii]   = 2*ui*vi - 4 -2*gamma - a
+        ii += 1
+        # dv'(i)/du(i)
+        dfdy.nzval[ii] = 3 - 2*ui*vi
+        ii +=1
+        if  iu+2<=2N
+            # du'(i+1)/du(i)
+            dfdy.nzval[ii] = gamma
+            ii +=1
+        end
+        
+        # iv: column belonging to d/dv(i)
+        if iv-2>=1
+            # dv'(i-1)/dv(i)
+            dfdy.nzval[ii] = gamma
+            ii +=1
+        end
+        # du'(i)/dv(i)
+        dfdy.nzval[ii] = ui^2
+        ii +=1
+        # dv'(i)/dv(i)
+        dfdy.nzval[ii]   = -ui^2 -2*gamma - a
+        ii +=1
+        if  iu+2<=2N
+            # du'(i+1)/dv(i)==0
+            dfdy.nzval[ii] = 0
+            ii +=1
+        end
+        if iv+2<=2N
+            # dv'(i+1)/dv(i)
+            dfdy.nzval[ii] = gamma
+            ii +=1
+        end
+    end
+    return nothing
+end
+function jac!(;T_::Type=T, dof_=dof)
+    B = (ones(T_,dof_-2), ones(T_,dof_-1), ones(T_,dof_), ones(T_,dof_-1), ones(T_,dof_-2))        
+    spdiagm(B, (-2,-1,0,1,2))
+end
+function jac(t, y, ydot, a)
+    J = jac!()
+    jac!(J, y, ydot, a)
+    return J
+end
+function jac_ex(t, y)
+    return jac(t, y, y, 0)
+end
+
+refsol = T[0.9949197002317599, 3.0213845767604077, 0.9594350193986054, 3.0585989778165419, 0.9243010095428502, 3.0952478919989637, 0.8897959106772672,
+           3.1310118289054731, 0.8561653620284367, 3.1656101198770159, 0.8236197147449046, 3.1988043370624344, 0.7923328094811884, 3.2303999530641514,
+           0.7624421042573115, 3.2602463873623941, 0.7340499750795348, 3.2882356529108807, 0.7072259700779899, 3.3142998590079271, 0.6820097782458483,
+           3.3384078449410937, 0.6584146743834650, 3.3605612157873943, 0.6364312187752559, 3.3807900316323134, 0.6160310186921587, 3.3991483695914764,
+           0.5971703941198909, 3.4157099395342736, 0.5797938277687891, 3.4305638938070224, 0.5638371159206763, 3.4438109320334580, 0.5492301695479158,
+           3.4555597666485198, 0.5358994429426996, 3.4659239846027008, 0.5237699892215797, 3.4750193162238476, 0.5127671585747183, 3.4829613034792271,
+           0.5028179665048467, 3.4898633463634923, 0.4938521662914935, 3.4958350971335204, 0.4858030633656755, 3.5009811668111510, 0.4786081100251151,
+           3.5054001059792705, 0.4722093177200750, 3.5091836216744015, 0.4665535216425440, 3.5124159935026285, 0.4615925290790646, 3.5151736544621075,
+           0.4572831793403656, 3.5175249049438184, 0.4535873393501199, 3.5195297317024448, 0.4504718553589467, 3.5212397070273984, 0.4479084778719241,
+           3.5226979467564341, 0.4458737738041973, 3.5239391090719634, 0.4443490371324889, 3.5249894191569453, 0.4433202068820853, 3.5258667077466495,
+           0.4427777991494095, 3.5265804544017270, 0.4427168579654424, 3.5271318289682063, 0.4431369281018266, 3.5275137272135266, 0.4440420513508381,
+           3.5277107990730161, 0.4454407863109616, 3.5276994703501980, 0.4473462502188303, 3.5274479611304068, 0.4497761798232572, 3.5269163066324394,
+           0.4527530066369863, 3.5260563887768472, 0.4563039400688689, 3.5248119894251024, 0.4604610498812091, 3.5231188790654930, 0.4652613370907894,
+           3.5209049576992761, 0.4707467798082714, 3.5180904678044698, 0.4769643375804777, 3.5145883024867057, 0.4839658945842979, 3.5103044351908528,
+           0.4918081185812277, 3.5051385005173827, 0.5005522089940899, 3.4989845585737802, 0.5102635039989190, 3.4917320776245013, 0.5210109134090777,
+           3.4832671712209993, 0.5328661417420937, 3.4734741260299615, 0.5459026646938675, 3.4622372546582585, 0.5601944229089820, 3.4494431032230182,
+           0.5758142001453760, 3.4349830354873361, 0.5928316594749734, 3.4187562033108012, 0.6113110218368440, 3.4006728962523969, 0.6313083867734524,
+           3.3806582409098729, 0.6528687160104193, 3.3586561928427350, 0.6760225267555723, 3.3346337311179157, 0.7007823726569721, 3.3085851288057930,
+           0.7271392249346637, 3.2805361342349380, 0.7550589020044152, 3.2505478606008622, 0.7844787296769868, 3.2187201496972175, 0.8153046416214843,
+           3.1851941538893653, 0.8474089465959840, 3.1501538739882800, 0.8806289904192589, 3.1138264039027113, 0.9147669230929857, 3.0764806689389470,
+           0.9495907429372025, 3.0384245041548366, 0.9848367306701233] # reference sol from Hairer has 1.36 scd
+refsolinds = collect(1:7:dof) # refsol is only at components 1:7:2N
+
+
+ic = zeros(T,2N)
+ic[1:2:2N] = u0(x)
+ic[2:2:2N] = v0(x)
+tspan = T[0.0, 10.0] # integration interval
+tspan0 = T[0.0, 0.001] # integration interval
+
+#t,ydassl = dasslSolve(fn, ic, tspan, jacobian=jac)
+
+## ode23s
+# tout, yout = ode23s(fn_ex, ic, tspan0, jacobian=jac_ex)
+# @time tout, yout = ode23s(fn_ex, ic, tspan, jacobian=jac_ex)
+# @show norm(abs(yout[end][refsolinds]-refsol)./refsol, Inf)
+
+# ## DASSL
+# t,ydassl = dasslSolve(fn, ic, tspan0, jacobian=jac)
+# @time t,ydassl = dasslSolve(fn, ic, tspan, jacobian=jac)
+# @show norm(abs(ydassl[end][refsolinds]-refsol)./refsol, Inf)
+
+## ROSW
+
+# # No jac works but is slower and gives much higher precision than
+# # other methods.
+# tout, yout = ode_rosw(fn!, ic, tspan0)
+# @time tout, yout = ode_rosw(fn!, ic, tspan)
+# @show norm(abs(yout[end][refsolinds]-refsol)./refsol, Inf)
+
+# Analytic Jacobian
+tout, yout = ode_rosw(fn!, ic, tspan0, jacobian=(jac!, jac!()))
+@time tout, yout = ode_rosw(fn!, ic, tspan, jacobian=(jac!, jac!()))
+@show norm(abs(yout[end][refsolinds]-refsol)./refsol, Inf)
+
+# Numerical colored Jacobian
+tout, yout = ode_rosw(fn!, ic, tspan0, jacobian=jac!())
+@time tout, yout = ode_rosw(fn!, ic, tspan, jacobian=jac!())
+@show norm(abs(yout[end][refsolinds]-refsol)./refsol, Inf)
+
diff --git a/examples/hb1dae.jl b/examples/hb1dae.jl
index 461b7e40d..e6902685f 100644
--- a/examples/hb1dae.jl
+++ b/examples/hb1dae.jl
@@ -54,12 +54,12 @@ yr = hcat(yrosw...);
 println("Relative difference between DASSL vs ROSW, no Jac:")
 println(abs(ydassl[end]-yrosw[end])./abs(ydassl[end]))
 
-println("Using fixed step ROSW")
-tspan = t
-t,yrosw = ode_rosw_fixed(hb1dae!, y0, tspan)
-@time t,yrosw = ode_rosw_fixed(hb1dae!, y0, tspan)
-println("Relative difference between DASSL vs fixed step, no Jac:")
-println(abs(ydassl[end]-yrosw[end])./abs(ydassl[end]))
+# println("Using fixed step ROSW")
+# tspan = t
+# t,yrosw = ode_rosw_fixed(hb1dae!, y0, tspan)
+# @time t,yrosw = ode_rosw_fixed(hb1dae!, y0, tspan)
+# println("Relative difference between DASSL vs fixed step, no Jac:")
+# println(abs(ydassl[end]-yrosw[end])./abs(ydassl[end]))
 
 # with Jacobian
 tspan = [0, 4e6]
diff --git a/src/rosw.jl b/src/rosw.jl
index 04680c420..bd36d3c34 100644
--- a/src/rosw.jl
+++ b/src/rosw.jl
@@ -26,10 +26,12 @@ export ode_rosw, ode_rosw_fixed, ode_rodas3
 # The method used here uses transformed equations as done in PETSc.
 
 
-rosw_poststep(xtrial, err, t, dt, steps) = (showcompact(t), print(", ") ,
-                                            showcompact(dt), print(", "),
-                                            showcompact(err), print(", "),
-                                            println(steps))
+# rosw_poststep(xtrial, err, t, dt, steps) = (showcompact(t), print(", ") ,
+#                                             showcompact(dt), print(", "),
+#                                             showcompact(err), print(", "),
+#                                             println(steps))
+
+rosw_poststep(xtrial, err, t, dt, steps) = nothing
 
 # Tableaus
 ##########
@@ -132,13 +134,13 @@ const bt_ros_rodas3 = TableauRosW(:ros_rodas3, (3,4), Rational{Int},
 ###################
 # Fixed step solver
 ###################
-ode_rosw_fixed(fn, x0, tspan; kwargs...) = oderosw_fixed(fn, x0, tspan, bt_ros34pw2, kwargs...)
-function oderosw_fixed{N,S}(fn, x0::AbstractVector, tspan,
+ode_rosw_fixed(fn!, x0, tspan; kwargs...) = oderosw_fixed(fn!, x0, tspan, bt_ros34pw2, kwargs...)
+function oderosw_fixed{N,S}(fn!, x0::AbstractVector, tspan,
                             btab::TableauRosW{N,S};
-                            jacobian=numerical_jacobian(fn, 1e-6, 1e-6)
+                            jacobian=numerical_jacobian(fn!, 1e-6, 1e-6)
                             )
     # TODO: refactor with oderk_fixed
-    Et, Exf, Tx, btab = make_consistent_types(fn, x0, tspan, btab)
+    Et, Exf, Tx, btab = make_consistent_types(fn!, x0, tspan, btab)
     btab = tabletransform(btab)
     dof = length(x0)
 
@@ -155,26 +157,26 @@ function oderosw_fixed{N,S}(fn, x0::AbstractVector, tspan,
     u = zeros(Exf,dof)    # work vector 
     udot = zeros(Exf,dof) # work vector 
     
-    # allocate!(ks, x0, dof) # no need to allocate as fn is not in-place
+    # allocate!(ks, x0, dof) # no need to allocate as fn! is not in-place
     xtmp = similar(x0, Exf, dof)
     for i=1:length(tspan)-1
         dt = tspan[i+1]-tspan[i]
-        rosw_step!(xs[i+1], fn, jacobian, xs[i], dt, dof, btab,
+        rosw_step!(xs[i+1], fn!, jacobian, xs[i], dt, dof, btab,
                    k, jac_store, ks, u, udot)
     end
     return tspan, xs
 end
 
 
-ode_rosw(fn, x0, tspan;kwargs...) = oderosw_adapt(fn, x0, tspan, bt_ros34pw2; kwargs...)
-ode_rodas3(fn, x0, tspan;kwargs...) = oderosw_adapt(fn, x0, tspan, bt_ros_rodas3; kwargs...)
-function oderosw_adapt{N,S}(fn, x0::AbstractVector, tspan, btab::TableauRosW{N,S};
+ode_rosw(fn!, x0, tspan;kwargs...) = oderosw_adapt(fn!, x0, tspan, bt_ros34pw2; kwargs...)
+ode_rodas3(fn!, x0, tspan;kwargs...) = oderosw_adapt(fn!, x0, tspan, bt_ros_rodas3; kwargs...)
+function oderosw_adapt{N,S}(fn!, x0::AbstractVector, tspan, btab::TableauRosW{N,S};
                             reltol = 1.0e-5, abstol = 1.0e-8,
                             norm=Base.norm,
                             minstep=abs(tspan[end] - tspan[1])/1e18,
                             maxstep=abs(tspan[end] - tspan[1])/2.5,
                             initstep=0.,
-                            jacobian=numerical_jacobian(fn, reltol, abstol)
+                            jacobian=numerical_jacobian(fn!, reltol, abstol)
 #                            points=:all
                             )
     # TODO: refactor with oderk_adapt
@@ -188,8 +190,8 @@ function oderosw_adapt{N,S}(fn, x0::AbstractVector, tspan, btab::TableauRosW{N,S
         points=:all
     end
     ## Figure types
-    fn_expl = (t,x)->(out=similar(x); fn(out, x, x*0); out)
-    Et, Exf, Tx, btab = make_consistent_types(fn, x0, tspan, btab)
+    fn_expl = (t,x)->(out=similar(x); fn!(out, x, x*0); out)
+    Et, Exf, Tx, btab = make_consistent_types(fn!, x0, tspan, btab)
     btab = tabletransform(btab)
     # parameters
     order = minimum(btab.order)
@@ -211,10 +213,13 @@ function oderosw_adapt{N,S}(fn, x0::AbstractVector, tspan, btab::TableauRosW{N,S
     k = Array(Tx, S) # stage variables
     allocate!(k, x0, dof)
     ks = zeros(Exf,dof) # work vector for one k
-    if isa(jacobian, Tuple)
+    if isa(jacobian, Tuple) # Jacobian function and sparse storage provided
         jac_store = jacobian[2]
         jacobian = jacobian[1]
-    else
+    elseif isa(jacobian, SparseMatrixCSC) # Sparse storage provided, make colored Jacobian (TODO)
+        jac_store = jacobian
+        jacobian = numerical_jacobian(fn!, reltol, abstol) #, jac_store, jac_store)
+    else # nothing provided, make dense numerical jacobian
         jac_store = zeros(Exf,dof,dof) # Jacobian storage
     end
     u = zeros(Exf,dof)    # work vector 
@@ -247,11 +252,11 @@ function oderosw_adapt{N,S}(fn, x0::AbstractVector, tspan, btab::TableauRosW{N,S
     steps = [0,0]  # [accepted, rejected]
 
     ## Integration loop
-    laststep = false
+    islaststep = abs(t+dt-tend)<=eps(tend) ? true : false
     timeout = 0 # for step-control
     iter = 2 # the index into tspan and xs
     while true
-        rosw_step!(xtrial, fn, jacobian, x, dt, dof, btab,
+        rosw_step!(xtrial, fn!, jacobian, x, dt, dof, btab,
                    k, jac_store, ks, u, udot)
         # Completion again for embedded method, see line 927 of
         # http://www.mcs.anl.gov/petsc/petsc-current/src/ts/impls/rosw/rosw.c.html#TSROSW
@@ -276,13 +281,13 @@ function oderosw_adapt{N,S}(fn, x0::AbstractVector, tspan, btab::TableauRosW{N,S
             # f1 = fn_expl(t+dt, xtrial)
             if points==:specified
                 # interpolate onto given output points
-                while iter-1<nsteps_fixed && (tdir*tspan[iter]<tdir*(t+dt) || laststep) # output at all new times which are < t+dt
+                while iter-1<nsteps_fixed && (tdir*tspan[iter]<tdir*(t+dt) || islaststep) # output at all new times which are < t+dt
                     error("FIXME: not implemented yet")
                     hermite_interp!(xs[iter], tspan[iter], t, dt, x, xtrial, f0, f1) # TODO: 3rd order only!
                     iter += 1
                 end
             else
-                # first interpolate onto given output points
+                # First interpolate onto given output points
                 while iter_fixed-1<nsteps_fixed && tdir*t<tdir*tspan_fixed[iter_fixed]<tdir*(t+dt) # output at all new times which are < t+dt
                     error("FIXME: not implemented yet")
                     xout = hermite_interp(tspan_fixed[iter_fixed], t, dt, x, xtrial, f0, f1)
@@ -299,7 +304,7 @@ function oderosw_adapt{N,S}(fn, x0::AbstractVector, tspan, btab::TableauRosW{N,S
             # k[1] = f1 # load k[1]==f0 for next step
 
             # Break if this was the last step:
-            laststep && break
+            islaststep && break
 
             # Swap bindings of x and xtrial, avoids one copy
             x, xtrial = xtrial, x
@@ -311,14 +316,14 @@ function oderosw_adapt{N,S}(fn, x0::AbstractVector, tspan, btab::TableauRosW{N,S
             # Hit end point exactly if next step within 1% of end:
             if tdir*(t+dt*1.01) >= tdir*tend
                 dt = tend-t
-                laststep = true # next step is the last, if it succeeds
+                islaststep = true # next step is the last, if it succeeds
             end
         elseif abs(newdt)<minstep  # minimum step size reached, break
-            println(xerr)
+            println("max err=$(maximum(xerr)), mean err=$(mean(xerr))")
             println("Warning: dt < minstep.  Stopping.")
             break
         else # redo step with smaller dt
-            laststep = false
+            islaststep = false
             steps[2] +=1
             dt = newdt
             timeout = timeout_const
@@ -345,8 +350,9 @@ function rosw_step!{N,S}(xtrial, g!, gprime!, x, dt, dof, btab::TableauRosW_T{N,
         # first step of Newton iteration with guess ks==0
         #        k[s][:] = ks - jac\k[s] # in-place A_ldiv_B!(jac, k[s])
         # TODO: add option to do more Newton iterations
-        A_ldiv_B!(jac, k[s])  # TODO: see whether this is impacted by https://github.com/JuliaLang/julia/issues/10787
+        k[s] = A_ldiv_B!(jac, k[s])  # TODO: see whether this is impacted by https://github.com/JuliaLang/julia/issues/10787
                               # and https://github.com/JuliaLang/julia/issues/11325
+#        k[s][:]  = jac_store\k[s]
         for d=1:dof
             k[s][d] = ks[d]-k[s][d]
         end
@@ -356,7 +362,8 @@ function rosw_step!{N,S}(xtrial, g!, gprime!, x, dt, dof, btab::TableauRosW_T{N,
     # last step:
     #    k[S][:] = ks - jac\k[S]
     # TODO: add option to do more Newton iterations
-    A_ldiv_B!(jac, k[S])
+    k[S] = A_ldiv_B!(jac, k[S])
+#    k[S][:]  = jac_store\k[S]    
     for d=1:dof
         k[S][d] = ks[d]-k[S][d]
     end
@@ -408,7 +415,13 @@ function hprime!(res, xi, gprime!, u, udot, btab, dt)
     # The res will hold part of the LU factorization.  However use the
     # returned LU-type.
     gprime!(res, u, udot, 1./(dt*btab.γii))
-    return lufact!(res)
+    if true # issparse(res)
+        # For issues with sparse https://github.com/JuliaLang/julia/issues/7488
+        # Basically lufact! destroys res and it cannot be used anymore later.
+        return lufact(res)
+    else
+        return lufact!(res)
+    end
 end
 
 # Copied & adapted from DASSL:
@@ -416,12 +429,11 @@ end
 # finite differences
 function numerical_jacobian(fn!, reltol, abstol)
     function numjac!(jac, y, dy, a)
-        # TODO use coloring
-        
+        # Numerical Jacobian, finite differences, no coloring.
         ep      = eps(1e6)   # this is the machine epsilon # TODO: better value?
         h       = 1/a           # h ~ 1/a
         wt      = reltol*abs(y).+abstol
-        # delta for approximation of jacobian.  I removed the
+        # Delta for approximation of Jacobian.  I removed the
         # sign(h_next*dy0) from the definition of delta because it was
         # causing trouble when dy0==0 (which happens for ord==1)
         edelta  = max(abs(y),abs(h*dy),wt)*sqrt(ep)
@@ -437,7 +449,7 @@ function numerical_jacobian(fn!, reltol, abstol)
             tmp1[i] += edelta[i]
             tmp2[i] += a*edelta[i]
             fn!(tmp, tmp1, tmp2)
-            if isa(jac, SparseMatrixCSC)
+            if false #isa(jac, SparseMatrixCSC)
                 for nz in nzrange(jac,i)
                     nonzeros(jac)[nz] = tmpy[rowvals(jac)[nz]]
                 end
@@ -468,9 +480,7 @@ function numerical_jacobian(fn!, reltol, abstol, JPattern1::SparseMatrixCSC, JPa
         ep      = eps(1e6)   # this is the machine epsilon # TODO: better value?
         h       = 1/a           # h ~ 1/a
         wt      = reltol*abs(y).+abstol
-        # delta for approximation of jacobian.  I removed the
-        # sign(h_next*dy0) from the definition of delta because it was
-        # causing trouble when dy0==0 (which happens for ord==1)
+
         edelta  = max(abs(y),abs(h*ydot),wt)*sqrt(ep)
 
         dof = length(y)
@@ -486,7 +496,7 @@ function numerical_jacobian(fn!, reltol, abstol, JPattern1::SparseMatrixCSC, JPa
             fn!(res1, dy, ydot)
             fn!(res2, y, dydot)
             for j=1:dof
-                res1[j] = ( res1[j]-res0[j] + a*(res2[j]-res0[j]) )/edelta[i]
+                res1[j] = ( res1[j]-res0[j] + a*(res2[j]-res0[j]) )/edelta[j]
             end
             recover_jac!(jac, res1, S, c) # recovers the columns of jac which correspond c
             remove_delta!(dy, S, c, edelta) # remove delta again