     if (p.beta!=0.) apar[0] = induct(y); //computes a_par and needs correct npe
 //     if (p.flrmode==1) invert_invgammaA(invgammaDIR,apar[1] ,apar[0] );             //chi= Gamma (Ni-1)
     if (p.flrmode==1) {
-        old_gammaApar.extrapolate( apar[1]);
-        std::vector<unsigned> number = multigrid.direct_solve( multi_invgammaDIR, apar[1], apar[0] , p.eps_gamma);
-        old_gammaApar.update( apar[1]);
-        if(  number[0] == invert_invgamma.get_max())
-            throw dg::Fail( p.eps_gamma); 
+        if (p.tau[1] == 0.) {
+          dg::blas1::axpby(1.0,apar[0],0.,apar[1]);
+        } 
+        else {
+            old_gammaApar.extrapolate( apar[1]);
+            std::vector<unsigned> number = multigrid.direct_solve( multi_invgammaDIR, apar[1], apar[0] , p.eps_gamma);
+            old_gammaApar.update( apar[1]);
+            if(  number[0] == invert_invgamma.get_max())
+                throw dg::Fail( p.eps_gamma); 
+        }
     }
     if (p.flrmode==0) dg::blas1::axpby(1.0,apar[0],0.,apar[1]);
 

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

    
    old_gammaN.extrapolate( chi);
    std::vector<unsigned> number = multigrid.direct_solve( multi_invgammaN, chi, y[1], p.eps_gamma); //chi= Gamma (Ni-1)
    old_gammaN.update( chi);
    if(  number[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]);
    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);   

template<class Geometry, class IMatrix, class Matrix, class container>
container& Asela<Geometry, IMatrix, Matrix,container>::compute_psi( container& potential)
{
    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]);        
