anchor.diff:      !init for diff
       call fluxinitcont(   numcells,phiDSAdiff,phistart ) !create phi=phistart
       iter=0
  -    do while (error>tol)
  +    flitermore = .true.
  +    do while (flitermore)
         iter=iter+1
         !Step 1, Sn solve for half step solution of flux
         call FESn(         numcells,phiDSAlold,phiDSArold,psil,psir,x,dx,qr,ql,sigt,siga,mu,&
  
positive.code:   real(8),allocatable :: phiHsl(:),phiHsr(:),phiresl(:),phiresr(:)
    real(8),allocatable :: phiDSAlold(:),phiDSArold(:),phiDSAdiff(:)
    real(8) :: alpha,beta
  
    !Read and Prepare parameters
    call input(            a,curRight,curLeft,constq,qtype,phistart,&
  
negative.code: 
        if(iter>1) call FESnerror( numcells,phiDSAl,phiDSAr,phiDSAlold,phiDSArold,error ) !Converge?
        call FESnNewToOld( numcells,phiDSAl,phiDSAr,phiDSAlold,phiDSArold )
        if(error<tol) flitermore = .false.
      enddo
  !    print *,"iteration: ",iter
  !    call cpu_time(toc)
  