 template<class Geometry, class IMatrix, class Matrix, class container>
 container& Asela<Geometry, IMatrix, Matrix,container>::compute_psi( container& potential)
 {
-    old_psi.extrapolate( phi[1]);
-    std::vector<unsigned> number = multigrid.direct_solve( multi_invgammaDIR, phi[1], potential, p.eps_gamma);
-    old_psi.update( phi[1]);
-    if(  number[0] == invert_invgamma.get_max())
-      throw dg::Fail( p.eps_gamma); 
+    if (p.tau[1] == 0.) {
+        dg::blas1::axpby( 1., potential, 0., phi[1]); 
+    } 
+    else {
+        old_psi.extrapolate( phi[1]);
+        std::vector<unsigned> number = multigrid.direct_solve( multi_invgammaDIR, phi[1], potential, p.eps_gamma);
+        old_psi.update( phi[1]);
+        if(  number[0] == invert_invgamma.get_max())
+        throw dg::Fail( p.eps_gamma); 
+    }
     poissonN.variationRHS(potential, omega); 
     dg::blas1::pointwiseDot(1.0, binv, binv, omega, 0.0, omega);        // omega = u_E^2 
     dg::blas1::axpby( 1., phi[1], -0.5, omega,phi[1]);        

template<class Geometry, class IMatrix, class Matrix, class container>
void Asela<Geometry, IMatrix, Matrix, container>::initializene( const container& src, container& target)
{ 
    std::vector<unsigned> number = multigrid.direct_solve( multi_invgammaN, target,src, p.eps_gamma);  //=ne-1 = Gamma (ni-1)  
    if(  number[0] == invert_invgamma.get_max())
      throw dg::Fail( p.eps_gamma);
}

template<class G, class IMatrix, class M, class V>

        multi_pol[u].set_chi( multi_chi[u]);
    }

    //Gamma N_i
    if (p.tau[1] == 0.) {
        dg::blas1::axpby( 1., y[1], 0.,chi); //chi = N_i - 1
    } 
    else {
        old_gammaN.extrapolate( chi);
        std::vector<unsigned> numberG = multigrid.direct_solve( multi_invgammaN, chi, y[1], p.eps_gamma); //chi= Gamma (Ni-1)
        old_gammaN.update( chi);
        if(  numberG[0] == invert_invgamma.get_max())
            throw dg::Fail( p.eps_gamma);
    }
    //rhs
    dg::blas1::axpby( -1., y[0], 1.,chi,chi);        //chi=  Gamma (n_i-1) - (n_e-1) = Gamma n_i - n_e
    //polarisation
    old_phi.extrapolate( phi[0]);
    std::vector<unsigned> number = multigrid.direct_solve( multi_pol, phi[0], chi, p.eps_pol);
    old_phi.update( phi[0]);
    if(  number[0] == invert_pol.get_max())
        throw dg::Fail( p.eps_pol);   
