     //compute chi and polarisation
     if( global) 
     {
-        exp( y, expy);
-        blas1::pointwiseDot( binv, expy[1], chi); //\chi = n_i
+        dg::blas1::transfer( y[1], chi);
+        dg::blas1::plus( chi, 1.); 
+        blas1::pointwiseDot( binv, chi, chi); //\chi = n_i
         blas1::pointwiseDot( binv, chi, chi); //\chi *= binv^2
-        //A = B;  
         if( global == 1)
             pol.set_chi( chi);
-        //compute omega
-        dg::blas1::transform( expy[0], expy[0], dg::PLUS<double>(-1)); //n_e -1
-        dg::blas1::transform( expy[1], omega, dg::PLUS<double>(-1)); //n_i -1
     }
     else
     {
         blas1::axpby( 1., y[1], 0., omega); //n_i = omega
     }
-    blas2::symv( w2d, omega, omega); 
 #ifdef DG_BENCHMARK
     t.toc();
     std::cout<< "Polarisation assembly took "<<t.diff()<<"s\n";
     t.tic();
 #endif 
-    //Attention!! gamma1 wants Dirichlet BC
-    unsigned number = pcg( gamma1, gamma_n, omega, v2d, eps_gamma);
+    invert_invgamma( gamma1, gamma_n, y[1]);
     if( number == pcg.get_max())
         throw Fail( eps_gamma);
 #ifdef DG_BENCHMARK


};

template< class M, class container>
ToeflR< M, container>::ToeflR( const dg::CartesianGrid2d& grid, double kappa, double nu, double tau, double eps_pol, double eps_gamma, int global ): 
    chi( grid.size(), 0.), omega(chi), gamma_n( chi), gamma_old( chi), 
    binv( evaluate( LinearX( kappa, 1.), grid)), 
    phi( 2, chi), phi_old( phi), dyphi( phi),
    expy( phi), dxy( expy), dyy( dxy), lapy( dyy),
    gamma1(  grid, -0.5*tau, dg::centered),
    arakawa( grid), 
    pol(     grid, not_normed, dg::centered), 
    pcg( omega, omega.size()), 
    w2d( create::weights(grid)), v2d( create::inv_weights(grid)), one( dg::evaluate( dg::one, grid)), 
    eps_pol(eps_pol), eps_gamma( eps_gamma), kappa(kappa), nu(nu), tau( tau), global( global),
    laplaceM( grid, normed, dg::centered)
    //A( grid, not_normed)
{
}

template< class M, class container>
const container& ToeflR<M, container>::compute_vesqr( const container& potential)
{
    assert( global);
    blas2::gemv( arakawa.dx(), potential, chi);
    blas2::gemv( arakawa.dy(), potential, omega);
    blas1::pointwiseDot( binv, chi, chi);
    blas1::pointwiseDot( binv, omega, omega);
    blas1::pointwiseDot( chi, chi, chi);
    blas1::pointwiseDot( omega, omega, omega);
    blas1::axpby( 1., chi, 1.,  omega);
    return omega;
}
template< class M, class container>
const container& ToeflR<M, container>::compute_psi( const container& potential)
{
    //compute Gamma phi[0]
    blas1::axpby( 2., phi[1], -1.,  phi_old[1]);
    phi[1].swap( phi_old[1]);

    blas2::symv( w2d, potential, omega);
#ifdef DG_BENCHMARK
    Timer t;
    t.tic();
#endif //DG_BENCHMARK
    unsigned number = pcg( gamma1, phi[1], omega, v2d, eps_gamma);
    if( number == pcg.get_max())
        throw Fail( eps_gamma);
#ifdef DG_BENCHMARK
    std::cout << "# of pcg iterations for psi \t"<< number << "\t";
    t.toc();
    std::cout<< "took \t"<<t.diff()<<"s\n";
#endif //DG_BENCHMARK
    //now add -0.5v_E^2
    if( global)
    {
        blas1::axpby( 1., phi[1], -0.5, compute_vesqr( potential), phi[1]);
    }
    return phi[1];
}


//computes and modifies expy!!
template<class M, class container>
const container& ToeflR< M, container>::polarisation( const std::vector<container>& y)
{
    //USE INVERT CLASS 
    //extrapolate phi and gamma_n
    blas1::axpby( 2., phi[0], -1.,  phi_old[0]);
    blas1::axpby( 2., gamma_n, -1., gamma_old);
    //blas1::axpby( 1., phi[1], 0.,  phi_old[1]);
    //blas1::axpby( 0., gamma_n, 0., gamma_old);
    gamma_n.swap( gamma_old);
    phi[0].swap( phi_old[0]);

#ifdef DG_BENCHMARK
    Timer t; 
    t.tic();

#endif 
    if( global)
    {
        blas1::axpby( -1., y[0], 1., gamma_n, omega); //omega = a_i\Gamma n_i - n_e
        if( global == 2)
            divide( omega, expy[1], omega);
        blas2::symv( w2d, omega, omega);
