         std::vector<double> px = detail::coefficients( xn, g.n()), 
                             py = detail::coefficients( yn, g.n()),
                             pz = detail::coefficients( zn, 1 );
+        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> pxyz( g.n()*g.n());
         for(unsigned k=0; k<g.n(); k++)
             for( unsigned j=0; j<g.n(); j++)
-                pxyz[k*g.n()+j]= pz[0]*py[k]*px[j];
+                pxyz[k*g.n()+j]= pz[0]*pyF[k]*pxF[j];
 
         //...these are the matrix coefficients with which to multiply 
         if (globalbcz == dg::DIR)

    //assert( z.size() == g.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
