     curvZ( dg::evaluate(solovev::CurvatureZ(gp), g)),
     gradlnB( dg::evaluate(solovev::GradLnB(gp) , g)),
     iris( dg::evaluate( solovev::Pupil( gp), g)),
-    source( dg::evaluate( dg::Gaussian( p.R_0, 0, p.b, p.b, p.amp_source, 0), g)),
-    damping( dg::evaluate( solovev::Damping(gp ), g)), 
+//     source( dg::evaluate( dg::Gaussian( gp.R_0, 0, p.b, p.b, p.amp_source, 0), g)),
+//     source( dg::evaluate( solovev::Gradient(gp), g)),
+     source( dg::evaluate(solovev::TanhSource(gp, p.amp_source), g)),
+    damping( dg::evaluate( solovev::TanhDamping(gp ), g)), 
     phi( 2, chi), curvphi( phi), dzphi(phi), expy(phi),  
     dzy( 4, chi), curvy(dzy),
     A (dg::create::laplacianM_perp( g, dg::not_normed, dg::symmetric)),
-    dz(solovev::Field(gp), g),
+    dz(solovev::Field(gp), g, gp.rk4eps),
     arakawa( g), 
     pol(     g), 
     invert_pol( omega, omega.size(), p.eps_pol), 

        gp(gp),
        w3d_( dg::create::w3d(g)), v3d_(dg::create::v3d(g)),
        temp( g.size()),
        pupil_( dg::evaluate( solovev::Pupil( gp), g))
    {
        LaplacianM_perp = dg::create::laplacianM_perp( g, dg::normed, dg::symmetric);
        //LaplacianM_para = dg::create::laplacianM_parallel( g, dg::PER);
    }
    void operator()( const std::vector<container>& x, std::vector<container>& y)
    {

        {
            dg::blas2::gemv( LaplacianM_perp, x[i], temp);
            dg::blas2::gemv( LaplacianM_perp, temp, y[i]);
            dg::blas1::axpby( -p.nu_perp, y[i], 0., y[i]); // - nu_perp lapl_RZ (lapl_RZ (lnN,U)) //factor MISSING!?!
//             dg::blas2::gemv( LaplacianM_para, x[i], temp);
//             dg::blas1::axpby(  -p.nu_parallel, temp, 1., y[i]);
        }
        //add parallel resistivity
        std::vector<container>  expy(2);
