diff --git a/src/operators/diffusion/CMakeLists.txt b/src/operators/diffusion/CMakeLists.txt
index d772eef4bb8a7962de17236a348295358511eeec..7ab4da55c65507e22a6ce23bc8efd1d551343f16 100755
--- a/src/operators/diffusion/CMakeLists.txt
+++ b/src/operators/diffusion/CMakeLists.txt
@@ -4,10 +4,10 @@ SET( headers tnlLinearDiffusion.h
              tnlLinearDiffusion_impl.h
              tnlExactLinearDiffusion.h             
              tnlExactLinearDiffusion_impl.h
-	     tnlNonLinearDiffusion.h
-             tnlNonLinearDiffusion_impl.h
-             tnlExactNonLinearDiffusion.h             
-             tnlExactNonLinearDiffusion_impl.h )
+	     tnlNonlinearDiffusion.h
+             tnlNonlinearDiffusion_impl.h
+             tnlExactNonlinearDiffusion.h             
+             tnlExactNonlinearDiffusion_impl.h )
 
 SET( CURRENT_DIR ${CMAKE_SOURCE_DIR}/src/operators/diffusion )
 
diff --git a/src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator_impl.h b/src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator_impl.h
index 095f1d437a146391374e6742987b233a011ae94e..a42a70e4724a7ced98f59be5a39796a6baaa3e6f 100644
--- a/src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator_impl.h
+++ b/src/operators/diffusion/nonlinear-diffusion-operators/tnlOneSideDiffNonlinearOperator_impl.h
@@ -88,7 +88,7 @@ updateLinearSystem( const RealType& time,
 {
    const RealType aCoef = - tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * mesh.getHxSquareInverse() / 
                        operatorQ.getValue(mesh, mesh.template getCellNextToCell< -1 >( index ), coordinates, u, time );
-   const RealType bCoef = 1 + tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * ( mesh.getHxSquareInverse() / 
+   const RealType bCoef = tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * ( mesh.getHxSquareInverse() / 
                           operatorQ.getValue(mesh, index, coordinates, u, time ) + mesh.getHxSquareInverse() / 
                           operatorQ.getValue(mesh, mesh.template getCellNextToCell< -1 >( index ), coordinates, u, time ) );
    const RealType cCoef = - tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * mesh.getHxSquareInverse() / 
@@ -186,7 +186,7 @@ updateLinearSystem( const RealType& time,
                        operatorQ.getValue(mesh, mesh.template getCellNextToCell< 0,-1 >( index ), coordinates, u, time );
    const RealType bCoef = - tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * mesh.getHxSquareInverse() / 
                        operatorQ.getValue(mesh, mesh.template getCellNextToCell< -1,0 >( index ), coordinates, u, time );
-   const RealType cCoef = 1 + tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * ( mesh.getHySquareInverse() / 
+   const RealType cCoef = tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * ( mesh.getHySquareInverse() / 
                        operatorQ.getValue(mesh, index, coordinates, u, time ) + mesh.getHySquareInverse() / 
                        operatorQ.getValue(mesh, mesh.template getCellNextToCell< 0,-1 >( index ), coordinates, u, time )
                        + mesh.getHxSquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time ) + 
@@ -296,7 +296,7 @@ updateLinearSystem( const RealType& time,
                        operatorQ.getValue(mesh, mesh.template getCellNextToCell< 0,-1,0 >( index ), coordinates, u, time );
    const RealType cCoef = - tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * mesh.getHxSquareInverse() / 
                        operatorQ.getValue(mesh, mesh.template getCellNextToCell< -1,0,0 >( index ), coordinates, u, time );
-   const RealType dCoef = 1 + tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * ( mesh.getHySquareInverse() / 
+   const RealType dCoef = tau * operatorQ.getValueStriped(mesh, index, coordinates, u, time ) * ( mesh.getHySquareInverse() / 
                        operatorQ.getValue(mesh, index, coordinates, u, time ) + mesh.getHySquareInverse() / 
                        operatorQ.getValue(mesh, mesh.template getCellNextToCell< 0,-1,0 >( index ), coordinates, u, time )
                        + mesh.getHxSquareInverse() / operatorQ.getValue(mesh, index, coordinates, u, time ) +