anchor.diff:      ! TODO
   
       ! Copy the result
  -    do i = 1, n_cell
  -       state_out(i)   = state_in(i)
  -       state_out(i)%S = S(i)
  -       state_out(i)%T = T(i)
  -    end do
     end subroutine do_timestep
     
   end module mod_integrator
  
positive.code:   
  contains
    
    subroutine do_timestep (state_in, state_out)
      use mod_constants
      use mod_timestep
      implicit none
      
      type (state), intent(in), dimension(:)        :: state_in
      type (state), intent(out), dimension(:)       :: state_out
  
      real(kind = x_precision)                      :: dt, dx2, dx = 0
      integer                                       :: i, info
  
negative.code:     real(kind = x_precision), dimension(n_cell-1) :: diag_low, diag_up
      
      ! Get the timestep
      call timestep(state, dt)
  
      dx2 = dx**2
      
      ! Create the diagonals and copy T and S
      diag     = 1/dt + 2/dx2 * state(:)%nu / state(:)%x
      diag_up  = -1/dx2 * state(2:n_cell)%nu / state(:)%x**2
      diag_low = -1/dx2 * state(1:n_cell-1)%nu / state(:)%x**2
  
      diag_up(1)           = -1/dx2 * state(2)%nu / state(1)%x**2
      diag(n_cell)         = state(n_cell)%nu
      diag_low(n_cell - 1) = state(n_cell-1)%nu
      
      ! Solving for S
      call dgtsv(n_cell, 1, diag_low, diag, diag_up, state_out(:)%S, n_cell, info)
      
      if (info /= 0) then
         print *, "Ooops, something bad happened!"
  