     const dg::DVec solution = dg::evaluate( deriNEU, g3d);
     const dg::DVec solutionT = dg::evaluate( deriNEUT, g3d);
     const dg::DVec solutiondzTdz = dg::evaluate( deriNEUT2, g3d);
+    const dg::DVec bR_ = dg::evaluate( bR, g3d);
+    const dg::DVec bZ_ = dg::evaluate( bZ, g3d);
+    const dg::DVec bPhi_ = dg::evaluate(bPhi, g3d);
     dz.set_boundaries( dg::PER, 0, 0);
 
     dz( function, derivative); //dz(f)

    private:
    double R_0, I_0;
};
//psi = cos(0.5*pi*(R-R_0))*cos(0.5*pi*Z)
struct Field
{

    private:
    double R_0, I_0;
};

double funcNEU(double R, double Z, double phi)
{
