         //evaluate 2d Legendre polynomials at (xn, yn)...
         std::vector<double> px = detail::coefficients( xn, g.n()), 
                             py = detail::coefficients( yn, g.n());
+        std::vector<double> pxF(g.n(),0), pyF(g.n(), 0);
+        for( unsigned l=0; l<g.n(); l++)
+            for( unsigned k=0; k<g.n(); k++)
+            {
+                pxF[l]+= px[k]*forward(k,l);
+                pyF[l]+= py[k]*forward(k,l);
+            }
         std::vector<double> pxy( g.n()*g.n());
         //these are the matrix coefficients with which to multiply 
-        for(unsigned k=0; k<py.size(); k++)
-            for( unsigned l=0; l<px.size(); l++)
-                pxy[k*px.size()+l]= py[k]*px[l];
+        for(unsigned k=0; k<pyF.size(); k++)
+            for( unsigned l=0; l<pxF.size(); l++)
+                pxy[k*px.size()+l]= pyF[k]*pxF[l];
         if (globalbcz == dg::DIR)
         {
             if ( x[i]==g.x0() || x[i]==g.x1()  || y[i]==g.y0()  || y[i]==g.y1())

        //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 
        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;
}
/**
 * @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());

    dg::Operator<double> forward( g.dlt().forward());
    int number = 0;
    for( unsigned i=0; i<x.size(); i++)
    {
