     if (p.flrmode == 1)
     {
         dg::blas1::axpby( p.beta/p.mu[0], npe[0], 0., chi); //chi = beta/mu_e N_e
-//         dg::blas1::axpby(- p.beta/p.mu[1],  npe[1], 1., chi); //chi =beta/mu_e N_e-beta/mu_i  N_i
         maxwell.set_chi(chi);
         dg::blas1::pointwiseDot( npe[1], y[3], chi);               //lambda = N_i w_i
         invert_invgammaNW(invgammaDIR,lambda,chi);             //chi= Gamma (Ni wi)
         dg::blas1::pointwiseDot( npe[0], y[2], chi);                 //chi     = n_e w_e
         dg::blas1::axpby( -1.,lambda , 1., chi);  //chi = - Gamma (n_i w_i) + n_e w_e
-        //maxwell = (lap_per - beta*(N_i/hatmu_i - n_e/hatmu_e)) A_parallel 
-        //chi=n_e w_e -N_i w_i
+        //maxwell = (lap_per + beta*( n_e/mu_e)) A_parallel 
+        //chi=n_e w_e -Gamma (N_i w_i )
         unsigned number = invert_maxwell( maxwell, apar[0], chi); //omega equals a_parallel
         if( number == invert_maxwell.get_max())
             throw dg::Fail( p.eps_maxwell);

    dg::blas1::transfer( dg::evaluate( dg::zero, g), lambda ); 
    dg::blas1::transfer( dg::evaluate( dg::one,  g), one); 
    phi.resize(2);apar.resize(2); phi[0] = phi[1] = apar[0]=apar[1] =  chi;
    curvapar.resize(2); curvapar[0] = curvapar[1] = chi;
    curvphi = curvkappaphi = npe = logn = u = u2 = un =  phi;
    dsphi = dslogn = dsun = dsu2 = poissonn = poissonu =  phi;
    poissonu2 =poissonw = poissonun = poissonlogn =poissonphi  = phi;

            vecdotnablaN(curvKappaX, curvKappaY, y[i], curvkappay[i]);           //K_kappa(N) = K_kappa(N-1)
            vecdotnablaDIR(curvKappaX, curvKappaY,  u[i], curvkappay[2+i]);      //K_kappa(U)
            vecdotnablaDIR(curvKappaX, curvKappaY, phi[i], curvkappaphi[i]);     //K_kappa(phi)
        
            if (p.bc==dg::DIR)
            {
                dg::blas1::pointwiseDot(u2[i],npe[i],chi); // N U^2
                vecdotnablaN(curvKappaX, curvKappaY, chi, lambda);      //K_kappa( N U^2)
                dg::blas1::pointwiseDot(u2[i],u[i],chi); // U^3
                vecdotnablaN(curvKappaX, curvKappaY, chi, omega);      //K_kappa( U^3)

                dg::blas1::axpby( -p.mu[i], lambda, 1., yp[i]);               //dtN = dtN - (hat(mu))  K_kappa(N U^2)
