     is a limiter and 0 if there isn't
     * @param field The field to integrate
     * @param grid The grid on which to operate
+    * @param hz 
     * @param eps Desired accuracy of runge kutta
     * @param limit Instance of the limiter class (Default is a limiter everywhere)
     * @param globalbcz Choose NEU or DIR. Defines BC in parallel on box

    container evaluate( BinaryOp f, UnaryOp g, unsigned p0, unsigned rounds);
    template< class BinaryOp, class UnaryOp>
    container evaluateAvg( BinaryOp f, UnaryOp g, unsigned p0, unsigned rounds);
    void eins(const Matrix& interp, const container& n, container& npe);
    void einsPlus( const container& n, container& npe);
    void einsMinus( const container& n, container& nme);
    void einsPlusT( const container& n, container& npe);
    void einsMinusT( const container& n, container& nme);
    void centeredT( const container& f, container& dzf);
    void forwardT( const container& f, container& dzf);
    void backwardT( const container& f, container& dzf);

    void symv( const container& f, container& dzTdzf);
    /**
     * @brief Returns the weights used to make the matrix symmetric 
     *

*
* @ingroup dz
* @tparam Matrix The matrix class of the interpolation matrix
* @tparam container The container-class on which the interpolation matrix operates on (does not need to be dg::HVec)
*/
template< class Matrix = dg::DMatrix, class container=thrust::device_vector<double> >
struct DZ
