diff --git a/src/DiffEqOperators.jl b/src/DiffEqOperators.jl
index 0af00efb2..9af8bd2a5 100644
--- a/src/DiffEqOperators.jl
+++ b/src/DiffEqOperators.jl
@@ -59,7 +59,7 @@ end
 export MatrixFreeOperator
 export AnalyticalJacVecOperator, JacVecOperator, getops
 export AbstractDerivativeOperator, DerivativeOperator,
-       CenteredDifference, UpwindDifference
+       CenteredDifference, UpwindDifference, NonLinearDiffusion
 export DirichletBC, Dirichlet0BC, NeumannBC, Neumann0BC, RobinBC, GeneralBC, MultiDimBC, PeriodicBC,
        MultiDimDirectionalBC, ComposedMultiDimBC
 export compose, decompose, perpsize
diff --git a/src/derivative_operators/derivative_operator.jl b/src/derivative_operators/derivative_operator.jl
index e87535959..fbfad008c 100644
--- a/src/derivative_operators/derivative_operator.jl
+++ b/src/derivative_operators/derivative_operator.jl
@@ -26,6 +26,32 @@ struct DerivativeOperator{T<:Real,N,Wind,T2,S1,S2<:SArray,T3,F} <: AbstractDeriv
     coeff_func              :: F
 end
 
+struct NonLinearDiffusion end
+
+function NonLinearDiffusion(second_differential_order::Int, first_differential_order::Int, approx_order::Int,
+                            p::AbstractVector{T}, q::AbstractVector{T}, dx::Union{T,AbstractVector{T}, Real},
+                            nknots::Int) where {T<:Real,N}
+    #p is given by bc1*k , k being the diffusion coefficient
+    #q is given by bc2*u , u being the desired function
+    #first differential is the inner one, second is the outer differential
+    @assert approx_order>1 "approximation_order must be greater than 1."
+    discretization = zeros(nknots)
+    
+    if first_differential_order > 0 
+        discretization += (CenteredDifference(first_differential_order,approx_order,dx,nknots)*q).*(CenteredDifference(second_differential_order,approx_order,dx,nknots)*p)
+    else 
+        discretization += q[2:(nknots + 1)].*(CenteredDifference(second_differential_order,approx_order,dx,nknots)*p) 
+    end
+    
+    for l = 1:(second_differential_order - 1)
+        discretization += binomial(second_differential_order,l)*(CenteredDifference(l + first_differential_order,approx_order,dx,nknots)*q).*(CenteredDifference(second_differential_order - l,approx_order,dx,nknots)*p)
+    end
+    
+    discretization += (CenteredDifference(first_differential_order + second_differential_order,approx_order,dx,nknots)*q).*p[2:(nknots + 1)]
+    return discretization
+    
+end
+
 struct CenteredDifference{N} end
 
 function CenteredDifference{N}(derivative_order::Int,