         //evaluate 2d Legendre polynomials at (xn, yn)...
         std::vector<double> px = detail::coefficients( xn, g.n());
         //...these are the matrix coefficients with which to multiply 
+        std::vector<double> pxF(px.size(),0);
+        for( unsigned l=0; l<g.n(); l++)
+            for( unsigned k=0; k<g.n(); k++)
+                pxF[l]+= px[k]*forward(k,l);
         unsigned col_begin = n*g.n();
-        detail::add_line( A, number, i,  col_begin, px);
-        //choose layout from comments
-    }
-    typedef cusp::coo_matrix<int, double, cusp::host_memory> Matrix;
-    
-    dg::Operator<double> forward( g.dlt().forward());
-    Matrix transformX = dg::tensor( g.N(), forward);
-
-    Matrix B( A);
-    cusp::multiply( A, transformX, B);
-    B.sort_by_row_and_column();
-    return B;
+        detail::add_line( A, number, i,  col_begin, pxF);
+    }
+    return A;
 }
 /**
  * @brief Create interpolation matrix

    assert( x.size() == y.size());
    cusp::coo_matrix<int, double, cusp::host_memory> A( x.size(), g.size(), x.size()*g.n()*g.n());

    int number = 0;
    for( unsigned i=0; i<x.size(); i++)
    {

        detail::add_line( A, number, i,  col_begin, g.n(), g.Nx(), pxyz);
        //choose layout from comments
    }
    return A;
}
/**
 * @brief Create interpolation between two grids
