// Newton pour l'oscillateur k - epsilon :: resolution en 1d
// - k'' + eps = 0
// - eps'' + alpha eps^2 / k = 0
// N = nb meshes
N = 4 ;
h = 1/N ;
alpha = 1.0 ;
EPS      = 1.0 ;
EPS_CONV = 1.0E-6 ;

// init vectors
XX     = ones(1,2*N) ;
XX_old = ones(1,2*N) ;
RHS    = ones(1,2*N) ;

// init jacobian
JJ     = zeros(2*N,2*N) ;

function [ent]=block_kk(i,j)
  ii = 0 ;
  jj = 0 ;
  ent = JJ(ii*N+i,jj*N+j) ;
endfunction
function [ent]=block_ke(i,j)
  ii = 0 ;
  jj = 1 ;
  ent = JJ(ii*N+i,jj*N+j) ;
endfunction
function [ent]=block_ek(i,j)
  ii = 1 ;
  jj = 0 ;
  ent = JJ(ii*N+i,jj*N+j) ;
endfunction
function [ent]=block_ee(i,j)
  ii = 1 ;
  jj = 1 ;
  ent = JJ(ii*N+i,jj*N+j) ;
endfunction

// diffusion operator
function [JJ]=assemble_diffusion(ii,jj)
  for i=1:N
    JJ(ii*N+i,jj*N+i) = JJ(ii*N+i,jj*N+i) + 2./h ;
  end
  for i=1:N-1
    JJ(ii*N+i,jj*N+i+1) = JJ(ii*N+i,jj*N+i+1) - 1./h ;
    JJ(ii*N+i+1,jj*N+i) = JJ(ii*N+i+1,jj*N+i) - 1./h ;
  end
endfunction

// reaction term
function [JJ]=assemble_reaction(r,ii,jj)
  for i=1:N
    JJ(ii*N+i,jj*N+i) = r * h ;
  end
endfunction

// --- assemble constant jacobian ---------------------------------------------
function []=assemble_constant_jacobian()
  // k-k : - Delta k
  assemble_diffusion(0,0) ;
  
  // k-epsilon : + epsilon
  assemble_reaction(1.0,0,1) ;
  
  // epsilon-epsilon : - Delta epsilon
  assemble_diffusion(1,1) ;
endfunction

// --- assemble nonlinear part of jacobian ------------------------------------
function []=assemble_nonlinear_jacobian()
  // epsilon-k : + alpha * epsilon^2 / k => - alpha * epsilon_old^2 / k_old^2 // OK
  for i=1:N
    block_ek(i,i) = - h * XX_old(N+i)^2 / XX_old(i)^2 ;
  end
  // epsilon-epsilon : + alpha * epsilon^2 / k => 2 * alpha * epsilon_old / k_old // OK
  for i=1:N
    block_ee(i,i) = 2.0 * h * XX_old(N+i) / XX_old(i) ;
  end
endfunction

// --- assemble source term ---------------------------------------------------
function []=assemble_rhs(xx)
  for i=1:N
   for j=1:N
    // - Delta k_old + epsilon_old
    xx(i)   = xx(i) + block_kk(i,j) * XX_old(j) ;
    xx(i)   = xx(i) + block_ke(i,j) * XX_old(N+j) ;
    // - Delta epsilon_old + alpha * epsilon^2 / k
    xx(N+i) = xx(N+i) + block_ek(i,j) * XX_old(j) ;
    xx(N+i) = xx(N+i) + block_ee(i,j) * XX_old(N+j) ;
    end
  end
endfunction

assemble_constant_jacobian() ;
while( EPS > EPS_CONV )

  RHS = zeros(1,2*N) ;

  // --- assemble jacobian ---
  assemble_nonlinear_jacobian()
    
  // --- assemble_source_term ---
  assemble_rhs(RHS)
  
  // --- solve ---
  XX = gmres(JJ,RHS') ;
  
  // --- compute residual --
  
  
  // --- save solution ---
  XX_old = XX ;
  
  XX
  
  // --- compute convergence ---
  
  EPS = 1.0E-7 ;

end