diff --git a/examples/inviscid-flow/1d/LaxFridrichs_impl.h b/examples/inviscid-flow/1d/LaxFridrichs_impl.h
index c1fa87b34c637da6edfb77e0fafda1a4c47c8399..6429a0c123c28b885c7955bf7f7f12d4110cd2df 100644
--- a/examples/inviscid-flow/1d/LaxFridrichs_impl.h
+++ b/examples/inviscid-flow/1d/LaxFridrichs_impl.h
@@ -40,6 +40,57 @@ operator()( const MeshFunction& u,
     static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); 
     static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); 
     const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); 
+      /*
+   //rho
+   const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); 
+   const IndexType& center = entity.getIndex(); 
+   const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); 
+   const IndexType& west = neighbourEntities.template getEntityIndex< -1 >();
+   typedef typename MeshType::Cell Cell;
+   const int size = mesh.template getEntitiesCount< Cell >()/5; 
+   return (0.5 * this->tau)( u[ west ] - 2.0 * u[ center ]  + u[ east ] ) 
+          - 0.5 * hxInverse * ( u[ west ] * u[ west + 3*size ] - u[ east ] * u[ east + 3*size ] );
+   */
+   /*
+   //rhoVel
+   const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); 
+   const IndexType& center = entity.getIndex(); 
+   const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); 
+   const IndexType& west = neighbourEntities.template getEntityIndex< -1 >();
+   typedef typename MeshType::Cell Cell;
+   const int size = mesh.template getEntitiesCount< Cell >()/5; 
+   return (0.5 * this->tau)( u[ west ] - 2.0 * u[ center ]  + u[ east ] ) 
+          - 0.5 * hxInverse * 
+          (( u[ west ] * u[ west + 2*size ] + u[ west + 3*size ] ) 
+          - (u[ east ] * u[ east + 2*size ] + u[ east + 3*size ] ));
+   */
+   /*
+   //energy
+   const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); 
+   const IndexType& center = entity.getIndex(); 
+   const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); 
+   const IndexType& west = neighbourEntities.template getEntityIndex< -1 >();
+   typedef typename MeshType::Cell Cell;
+   const int size = mesh.template getEntitiesCount< Cell >()/5; 
+   return (0.5 * this->tau)( u[ west ] - 2.0 * u[ center ]  + u[ east ] ) 
+          - 0.5 * hxInverse * 
+          (( u[ west ] + u[ west + 2*size ] ) * u[ west + size ]  
+          - (u[ east ] + u[ east + 2*size ] ) * u[ east + size ] );
+   */
+   /*
+   //vel
+   const IndexType& center = entity.getIndex(); 
+   typedef typename MeshType::Cell Cell;
+   const int size = mesh.template getEntitiesCount< Cell >()/5; 
+   return ( fu[ center - 2*size] - fu[ center - 3*size]);
+   */
+   /*
+   //pressure
+   const IndexType& center = entity.getIndex(); 
+   typedef typename MeshType::Cell Cell;
+   const int size = mesh.template getEntitiesCount< Cell >()/5; 
+   return ( this->gamma - 1 ) * ( fu[ center - 2*size ] - 0.5 * fu[ center - 3* size ] * fu[ center - 4*size]);
+   */
 
    const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -2 >(); 
    const IndexType& center = entity.getIndex(); 
diff --git a/examples/inviscid-flow/1d/euler.h b/examples/inviscid-flow/1d/euler.h
index cba3c12b65fa683360aa5076e861b7e69d6ec16a..1eeddb0179494f8382818fcc80579cf5ae4bfd4b 100644
--- a/examples/inviscid-flow/1d/euler.h
+++ b/examples/inviscid-flow/1d/euler.h
@@ -6,6 +6,10 @@
 #include <functions/tnlConstantFunction.h>
 #include "eulerProblem.h"
 #include "LaxFridrichs.h"
+#include "LaxFridrichsContinuity.h"
+#include "LaxFridrichsMomentum.h"
+#include "LaxFridrichsEnergy.h"
+
 #include "eulerRhs.h"
 #include "eulerBuildConfigTag.h"
 
diff --git a/examples/inviscid-flow/1d/eulerProblem.h b/examples/inviscid-flow/1d/eulerProblem.h
index 8c2e5020fc2417de5038ffbf73ba1f8f5b62c09a..0b9595b9bddb52736aeef3df88435052f3ec4cc1 100644
--- a/examples/inviscid-flow/1d/eulerProblem.h
+++ b/examples/inviscid-flow/1d/eulerProblem.h
@@ -25,6 +25,16 @@ class eulerProblem:
       using typename BaseType::MeshType;
       using typename BaseType::DofVectorType;
       using typename BaseType::MeshDependentDataType;
+      typedef tnlMeshFunction<Mesh,Mesh::getMeshDimensions(),RealType>MeshFunction;
+    //definition
+	tnlVector< RealType, DeviceType, IndexType > _uRho;
+	tnlVector< RealType, DeviceType, IndexType > _uRhoVelocity;
+	tnlVector< RealType, DeviceType, IndexType > _uEnergy;
+
+	tnlVector< RealType, DeviceType, IndexType > _fuRho;
+	tnlVector< RealType, DeviceType, IndexType > _fuRhoVelocity;
+	tnlVector< RealType, DeviceType, IndexType > _fuEnergy;
+
 
       static tnlString getTypeStatic();
       tnlSharedVector< RealType, DeviceType, IndexType > rho;
diff --git a/examples/inviscid-flow/1d/eulerProblem_impl.h b/examples/inviscid-flow/1d/eulerProblem_impl.h
index 06c647f9594586c9181c37b955ffe8d360456504..f4635838a3ca6a5cfd7f31ee17cb1c9ea5ad9d8f 100644
--- a/examples/inviscid-flow/1d/eulerProblem_impl.h
+++ b/examples/inviscid-flow/1d/eulerProblem_impl.h
@@ -106,7 +106,7 @@ setInitialCondition( const tnlParameterContainer& parameters,
    double eR = ( preR / (gamma - 1) ) + 0.5 * rhoR * velR * velR;
    double x0 = parameters.getParameter< double >( "riemann-border" );
    cout << gamma << " " << rhoL << " " << velL << " " << preL << " " << eL << " " << rhoR << " " << velR << " " << preR << " " << eR << " " << x0 << " " << gamma << endl;
-   int count = mesh.template getEntitiesCount< Cell >();
+   int count = mesh.template getEntitiesCount< Cell >()/3;
    this->rho.bind(dofs,0,count);
    this->rhoVel.bind(dofs,count,count);
    this->energy.bind(dofs,2 * count,count);
@@ -248,55 +248,65 @@ getExplicitRHS( const RealType& time,
     p this->pressure[]
     */
     typedef typename MeshType::Cell Cell;
-    int count = mesh.template getEntitiesCount< Cell >();
-    const RealType& size = 1;//mesh.template getSpaceStepsProducts< -1, 0 >();
-    for (long int i = 1; i < count - 1; i++)
-       _fu[i] = 1.0 / (2.0*tau) * (this->rho[i-1]+this->rho[i+1]-2.0*this->rho[i])-(1.0/(2.0 * size)) * (this->rho[i+1]
-       * this->velocity[i+1] - this->rho[i-1] * this->velocity[i-1]);
-       _fu[count-1] = _fu[count-2];
-    for (long int i = 1; i < count - 1; i++)
-       _fu[count + i] =
-       1.0 / (2.0 * tau) * (this->rhoVel[i+1] + this->rhoVel[i+1] - 
-       2.0 * this->rhoVel[i])-(1.0 / (2.0 * size)) * ((this->rhoVel[i+1]
-       * this->velocity[i + 1] + this->pressure[i + 1]) - (this->rhoVel[i-1] * this->velocity[i - 1] 
-       + this->pressure[i - 1]));
-       _fu[2*count-1] = _fu[2*count-2];
-    for (long int i = 1; i < count - 1; i++)
-       _fu[i + 2 * count] = 1.0 / (2.0*tau) * (this->energy[i-1]
-       + this->energy[i+1]-2.0*this->energy[i])
-       - (1.0/(2.0 * size)) * ((this->energy[i+1]
-       + this->pressure[i + 1]) * this->velocity[i + 1] - (this->energy[i-1] + this->pressure[i -1]) * this->velocity[i - 1]);
-       _fu[3 * count-1] = _fu[3 * count-2];
-    for (long int i = 0; i <= count; i++) //pressure
-       this->pressure[i] = (this->gamma - 1 ) * ( this->energy[i] - 0.5 * this->rhoVel[i] * this->velocity[i]);
-    for (long int i = 0; i <= count ; i++) //velocity
-       this->velocity[i] = this->rhoVel[i]/this->rho[i];
-     
-   /****
-    * If you use an explicit solver like tnlEulerSolver or tnlMersonSolver, you
-    * need to implement this method. Compute the right-hand side of
-    *
-    *   d/dt u(x) = fu( x, u )
-    *
-    * You may use supporting mesh dependent data if you need.
-    
+    int count = mesh.template getEntitiesCount< Cell >()/3;
+	//bind _u
+    this->_uRho.bind(_u,0,count);
+    this->_uRhoVelocity.bind(_u,count,count);
+    this->_uEnergy.bind(_u,2 * count,count);
+		
+	//bind _fu
+    this->_fuRho.bind(_u,0,count);
+    this->_fuRhoVelocity.bind(_u,count,count);
+    this->_fuEnergy.bind(_u,2 * count,count);
 
+   MeshFunctionType velocity( mesh, this->velocity );
+   MeshFunctionType pressure( mesh, this->pressure );
+	//rho
    this->bindDofs( mesh, _u );
    tnlExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;
-   MeshFunctionType u( mesh, _u ); 
-   MeshFunctionType fu( mesh, _fu ); 
+   MeshFunctionType uRho( mesh, _uRho ); 
+   MeshFunctionType fuRho( mesh, _fuRho );
+   diffrrentialOperatorRho.setTau(tau);
+   differentialOperatorRho.setVelocity(velocity) 
+   explicitUpdater.template update< typename Mesh::Cell >( time,
+                                                           mesh,
+                                                           this->differentialOperatorRho,
+                                                           this->boundaryCondition,
+                                                           this->rightHandSide,
+                                                           uRho,
+                                                           fuRho );
+   //rhoVelocity
+   MeshFunctionType uRhoVelocity( mesh, _uRhoVelocity ); 
+   MeshFunctionType fuRhoVelocity( mesh, _fuRhoVelocity );
+   diffrrentialOperatorRhoVelocity.setTau(tau);
+   differentialOperatorRhoVelocity.setVelocity(velocity)
+   differentialOperatorRhoVelocity.setPressure(pressure) 
+   explicitUpdater.template update< typename Mesh::Cell >( time,
+                                                           mesh,
+                                                           this->differentialOperatorRhoVelocity,
+                                                           this->boundaryCondition,
+                                                           this->rightHandSide,
+                                                           uRhoVelocity,
+                                                           fuRhoVelocity );
+   
+   //energy
+   MeshFunctionType uEnergy( mesh, _uEnergy ); 
+   MeshFunctionType fuEnergy( mesh, _fuEnergy );
+   diffrrentialOperatorEnergy.setTau(tau);
+   diffrrentialOperatorEnergy.setPressure(pressure);
+   diffrrentialOperatorEnergy.setVelocity(velocity); 
    explicitUpdater.template update< typename Mesh::Cell >( time,
                                                            mesh,
-                                                           this->differentialOperator,
+                                                           this->differentialOperatorEnergy,
                                                            this->boundaryCondition,
                                                            this->rightHandSide,
-                                                           u,
-                                                           fu );
+                                                           uEnergy,
+                                                           fuEnergy );
    tnlBoundaryConditionsSetter< MeshFunctionType, BoundaryCondition > boundaryConditionsSetter; 
    boundaryConditionsSetter.template apply< typename Mesh::Cell >( 
       this->boundaryCondition, 
       time + tau, 
-       u ); */
+       u );
  }
 
 template< typename Mesh,
diff --git a/examples/inviscid-flow/2d/LaxFridrichs_impl.h b/examples/inviscid-flow/2d/LaxFridrichs_impl.h
index 7f9c73731154d797334ccb86a26e7640743a7ae7..14c8a1132f782d694f15af8fea4214c57ea07a97 100644
--- a/examples/inviscid-flow/2d/LaxFridrichs_impl.h
+++ b/examples/inviscid-flow/2d/LaxFridrichs_impl.h
@@ -95,7 +95,7 @@ operator()( const MeshFunction& u,
 
 
 
-   const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); 
+   const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); 
    const IndexType& center = entity.getIndex(); 
    const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); 
    const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); 
diff --git a/examples/inviscid-flow/2d/euler.h b/examples/inviscid-flow/2d/euler.h
index ce35381216dbaaa170a3934247cb7e4d49a216d6..e5402806d7c71a4cd6efc9640683b78954562cc2 100644
--- a/examples/inviscid-flow/2d/euler.h
+++ b/examples/inviscid-flow/2d/euler.h
@@ -6,6 +6,10 @@
 #include <functions/tnlConstantFunction.h>
 #include "eulerProblem.h"
 #include "LaxFridrichs.h"
+#include "LaxFridrichsContinuity.h"
+#include "LaxFridrichsEnergy.h"
+#include "LaxFridrichsMomentumX.h"
+#include "LaxFridrichsMomentumY.h"
 #include "eulerRhs.h"
 #include "eulerBuildConfigTag.h"
 
diff --git a/examples/inviscid-flow/2d/eulerProblem.h b/examples/inviscid-flow/2d/eulerProblem.h
index 6a514c1da74cf0149feea852fdc227b08b1355cc..46f2b1f064eaa83a36d645f1e327d6a236a85994 100644
--- a/examples/inviscid-flow/2d/eulerProblem.h
+++ b/examples/inviscid-flow/2d/eulerProblem.h
@@ -25,6 +25,18 @@ class eulerProblem:
       using typename BaseType::MeshType;
       using typename BaseType::DofVectorType;
       using typename BaseType::MeshDependentDataType;
+
+    //definition
+	tnlVector< RealType, DeviceType, IndexType > _uRho;
+	tnlVector< RealType, DeviceType, IndexType > _uRhoVelocityX;
+	tnlVector< RealType, DeviceType, IndexType > _uRhoVelocityY;
+	tnlVector< RealType, DeviceType, IndexType > _uEnergy;
+
+	tnlVector< RealType, DeviceType, IndexType > _fuRho;
+	tnlVector< RealType, DeviceType, IndexType > _fuRhoVelocityX;
+	tnlVector< RealType, DeviceType, IndexType > _fuRhoVelocityY;
+	tnlVector< RealType, DeviceType, IndexType > _fuEnergy;
+
       tnlSharedVector< RealType, DeviceType, IndexType > rho;
       tnlSharedVector< RealType, DeviceType, IndexType > rhoVelX;
       tnlSharedVector< RealType, DeviceType, IndexType > rhoVelY;
diff --git a/examples/inviscid-flow/2d/eulerProblem_impl.h b/examples/inviscid-flow/2d/eulerProblem_impl.h
index aedadddb10efa559e2ec138974cd36f2b9368833..5cacac8daddca2b776675caf14bcfbe16fc7474f 100644
--- a/examples/inviscid-flow/2d/eulerProblem_impl.h
+++ b/examples/inviscid-flow/2d/eulerProblem_impl.h
@@ -229,91 +229,89 @@ getExplicitRHS( const RealType& time,
                 MeshDependentDataType& meshDependentData )
 {
     typedef typename MeshType::Cell Cell;
-    const RealType& cellSize = 1;// mesh.template getSpaceStepsProducts< -1, 0 >();
-    int size = mesh.template getEntitiesCount< Cell >()/4; 
-// prepsat na SWEN
-
-    for (long int j = 1; j < size - 1; j++)
-       {
-          for (long int i = 1; i < size - 1; i++)
-             {
-                _fu[j*size+i] = 1.0 / (4.0*tau) * 
-                (this->rho[j*size+i-1]+this->rho[j*size+i+1]+this->rho[(j-1)*size+i]+this->rho[(j+1)*size+i]-4.0*this->rho[j*size+i])
-                -(1.0/(2.0*cellSize))*(this->rho[j*size+i+1]*this->velocityX[j*size+i+1]-this->rho[j*size+i-1]*this->velocityX[j*size+i-1])
-                -(1.0/(2.0*cellSize))*(this->rho[(j+1)*size+i]*this->velocityY[(j+1)*size+i]-this->rho[(j-1)*size+i]*this->velocityY[(j-1)*size+i]);
-                _fu[pow(size,2)-size+i]=_fu[pow(size,2)-2*size+i];
-             };
-       _fu[j*size-1] = _fu[j*size-1];
-       };
-    for (long int j = 1; j < size - 1; j++)
-       {
-          for (long int i = 1; i < size - 1; i++)
-             {
-                _fu[pow(size,2) + j*size+i] = 1.0 / (4.0 * tau) *
-                (this->rhoVelX[j*size+i+1]+this->rhoVelX[j*size+i-1]+this->rhoVelX[(j-1)*size+i]+this->rhoVelX[(j+1)*size+i]-4.0*this->rhoVelX[j*size+i])
-               -(1.0/(2.0*cellSize))*((this->rhoVelX[j*size+i+1]*this->velocityX[j*size+i+1]+this->pressure[j*size+i+1])-(this->rhoVelX[j*size+i-1]*this->velocityX[j*size+i-1]+this->pressure[j*size+i-1]))
-               -(1.0/(2.0*cellSize))*((this->rhoVelX[(j+1)*size+i]*this->velocityY[(j+1)*size+i])-(this->rhoVelX[(j-1)*size+i]*this->velocityY[(j-1)*size+i]));
-               _fu[2*pow(size,2)-size+i]=_fu[2*pow(size,2)-2*size+i];
-             };
-       _fu[pow(size,2)+j*size] = _fu[pow(size,2)+j*size-1];
-       };
-    for (long int j = 1; j < size - 1; j++)
-       {
-          for (long int i = 1; i < size - 1; i++)
-             {
-                _fu[2*pow(size,2) + j*size+i] = 1.0 / (4.0 * tau) *
-                (this->rhoVelY[j*size+i+1]+this->rhoVelY[j*size+i-1]+this->rhoVelY[(j-1)*size+i]+this->rhoVelY[(j+1)*size+i]-4.0*this->rhoVelY[j*size+i])
-               -(1.0/(2.0*cellSize))*((this->rhoVelY[(j+1)*size+i]*this->velocityY[(j+1)*size+i]+this->pressure[(j+1)*size+i])-(this->rhoVelY[(j-1)*size+i]*this->velocityY[(j-1)*size+i]+this->pressure[(j-1)*size+i]))
-               -(1.0/(2.0*cellSize))*((this->rhoVelY[j*size+i+1]*this->velocityX[j*size+i+1])-(this->rhoVelY[j*size+i-1]*this->velocityX[j*size+i-1]));
-               _fu[3*pow(size,2)-size+i]=_fu[3*pow(size,2)-2*size+i];
-             };
-       _fu[2*pow(size,2)+j*size] = _fu[2*pow(size,2)+j*size-1];
-       };
-    for (long int j = 1; j < size - 1; j++)
-       {
-          for (long int i = 1; i < size - 1; i++)
-             {
-                _fu[3*pow(size,2) + j*size+i] = 1.0 / (4.0 * tau) *
-                (this->energy[j*size+i+1]+this->energy[j*size+i-1]+this->energy[(j-1)*size+i]+this->energy[(j+1)*size+i]-4.0*this->energy[j*size+i])
-       		-(1.0/(2.0*cellSize))*((this->energy[j*size+i+1]+this->pressure[j*size+i+1])*this->velocityX[j*size+i+1]-(this->energy[j*size+i-1]+this->pressure[j*size+i-1])*this->velocityX[j*size+i-1])
-        	-(1.0/(2.0*cellSize))*((this->energy[(j+1)*size+i]+this->pressure[(j+1)*size+i])*this->velocityY[(j+1)*size+i]-(this->energy[(j-1)*size+i]+this->pressure[(j-1)*size+i])*this->velocityY[(j-1)*size+i]);
-       		_fu[4*pow(size,2)-size+i] = _fu[4*pow(size,2)-2*size+i];
-             };
-       _fu[3*pow(size,2)+j*size] = _fu[3*pow(size,2)+j*size-1];
-       };
-
-    for (long int j = 1; j < size; j++) //pressure
-       for (long int i = 1; i < size; i++)
-          this->pressure[j*size+i] = (this->gamma - 1 ) * ( this->energy[j*size+i] - 0.5 * this->rho[j*size+i] * (pow(this->velocityX[j*size+i],2) + pow(this->velocityY[j*size+i],2)));
-    for (long int j = 1; j < size; j++) //velocityX
-       for (long int i = 1; i < size; i++)
-       this->velocityX[j*size+i] = this->rhoVelX[j*size+i]/this->rho[j*size+i];
-    for (long int j = 1; j < size; j++) //velocityY
-       for (long int i = 1; i < size; i++)
-       this->velocityY[j*size+i] = this->rhoVelY[j*size+i]/this->rho[j*size+i];
-    for (long int j = 1; j < size; j++) //velocity
-       for (long int i = 1; i < size; i++)
-       this->velocity[j*size+i] =sqrt(pow(velocityX[j*size+i],2)+pow(velocityY[j*size+i],2));
-
-/*
+    int count = mesh.template getEntitiesCount< Cell >()/4;
+	//bind _u
+    this->_uRho.bind(_u,0,count);
+    this->_uRhoVelocityX.bind(_u,count,count);
+    this->_uRhoVelocityY.bind(_u,2 * count,count);
+    this->_uEnergy.bind(_u,3 * count,count);
+		
+	//bind _fu
+    this->_fuRho.bind(_u,0,count);
+    this->_fuRhoVelocityX.bind(_u,count,count);
+    this->_fuRhoVelocityY.bind(_u,2 * count,count);
+    this->_fuEnergy.bind(_u,3 * count,count);
 
+   MeshFunctionType velocity( mesh, this->velocity );
+   MeshFunctionType velocityX( mesh, this->velocityX );
+   MeshFunctionType velocityY( mesh, this->velocityY );
+   MeshFunctionType pressure( mesh, this->pressure );
+   //rho
    this->bindDofs( mesh, _u );
    tnlExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;
-   MeshFunctionType u( mesh, _u ); 
-   MeshFunctionType fu( mesh, _fu ); 
+   MeshFunctionType uRho( mesh, _uRho ); 
+   MeshFunctionType fuRho( mesh, _fuRho );
+   diffrrentialOperatorRho.setTau(tau);
+   differentialOperatorRho.setVelocityX(velocityX);
+   differentialOperatorRho.setVelocityY(velocityY); 
    explicitUpdater.template update< typename Mesh::Cell >( time,
                                                            mesh,
-                                                           this->differentialOperator,
+                                                           this->differentialOperatorRho,
                                                            this->boundaryCondition,
                                                            this->rightHandSide,
-                                                           u,
-                                                           fu );
+                                                           uRho,
+                                                           fuRho );
+   //rhoVelocityX
+   MeshFunctionType uRhoVelocityX( mesh, _uRhoVelocityX ); 
+   MeshFunctionType fuRhoVelocityX( mesh, _fuRhoVelocityX );
+   diffrrentialOperatorRhoVelocityX.setTau(tau);
+   differentialOperatorRhoVelocityX.setVelocityX(velocityX);
+   differentialOperatorRhoVelocityX.setVelocityY(velocityY);
+   differentialOperatorRhoVelocityX.setPressure(pressure); 
+   explicitUpdater.template update< typename Mesh::Cell >( time,
+                                                           mesh,
+                                                           this->differentialOperatorRhoVelocityX,
+                                                           this->boundaryCondition,
+                                                           this->rightHandSide,
+                                                           uRhoVelocityX,
+                                                           fuRhoVelocityX );
+   //rhoVelocityY
+   MeshFunctionType uRhoVelocityY( mesh, _uRhoVelocityY ); 
+   MeshFunctionType fuRhoVelocityY( mesh, _fuRhoVelocityY );
+   diffrrentialOperatorRhoVelocityY.setTau(tau);
+   differentialOperatorRhoVelocityY.setVelocityX(velocityX);
+   differentialOperatorRhoVelocityY.setVelocityY(velocityY);
+   differentialOperatorRhoVelocityY.setPressure(pressure); 
+   explicitUpdater.template update< typename Mesh::Cell >( time,
+                                                           mesh,
+                                                           this->differentialOperatorRhoVelocityY,
+                                                           this->boundaryCondition,
+                                                           this->rightHandSide,
+                                                           uRhoVelocityY,
+                                                           fuRhoVelocityY );
+   
+   
+   //energy
+   MeshFunctionType uEnergy( mesh, _uEnergy ); 
+   MeshFunctionType fuEnergy( mesh, _fuEnergy );
+   diffrrentialOperatorEnergy.setTau(tau);
+   diffrrentialOperatorEnergy.setVelocityX(velocityX); 
+   diffrrentialOperatorEnergy.setVelocityY(velocityY); 
+   diffrrentialOperatorEnergy.setPressure(pressure);
+   explicitUpdater.template update< typename Mesh::Cell >( time,
+                                                           mesh,
+                                                           this->differentialOperatorEnergy,
+                                                           this->boundaryCondition,
+                                                           this->rightHandSide,
+                                                           uEnergy,
+                                                           fuEnergy );
+
+
    tnlBoundaryConditionsSetter< MeshFunctionType, BoundaryCondition > boundaryConditionsSetter; 
    boundaryConditionsSetter.template apply< typename Mesh::Cell >( 
       this->boundaryCondition, 
       time + tau, 
-       u ); 
-*/
+       u );
  }
 
 template< typename Mesh,