!
!  This code integrates a one-parameter flow on the line displaying 
!  two saddle-point bifurcations
!
!    dx      
!    -- = 1 + r*x + x**2     
!    dt    
!
!  The system displays two saddle-point bifurcations at 
!  r_c = -2 ( x* = 1 ) and r_c = 2 ( x* = -1 ) 
! 
!  For large r one has the following asymptotics :
!
!  x*+ ----> 1/r
!  x*- ----> -r
!
!  Instructions
!
!  - compile the code to produce an executable (.gfx extension)
!    
!          gfortran code.f90 -o code.gfx 
!  
!  - run the executable and redirect the standard output (i.e. 
!    screen on a file)
!
!          ./code.gfx > file_name.dat
!
!  - This will produce a two-column file (t,x), use gnuplot
!    to plot your results. Type gnuplot at the prompt in a 
!    terminal and then type at the gnuplot prompt
!
!          p 'file_name.dat' w l 
!
!    In gnuplot one uses abbreviations. For example
!
!    p   = plot
!    w l = with lines
!    pt 6 = pointtype 6 (open circles)
!    pt 7 = pointtype 7 (closed circles)
!
!    Refer to the website www.gnuplot.info for simple tutorials.
!
!================================================================
! Declarations
!
  integer, parameter :: kdp = kind(0.d0)
  real (kind=kdp) dt, T_tot, x, xt, ts, tk, L0, r
  integer N
!
! Parameters
!
  dt    = 0.001             ! Integrattion step
  dt2   = dt / 2.d0         ! half step
  T_tot = 10.               ! Total integration time (units of tau)
  N     = int ( T_tot/dt )  ! Number of integration steps
  ts    = 0.1               ! Sampling time (units of tau)
  
  tk    = 0.0               ! Initial value of running sampling time


  r     = 2.524367          ! The parameter of the one-dimensional flow
  x     = -20.              ! Initial value: look at your phase portrait
                            ! and bifurcation diagram to predict 
                            ! the asymptotic behaviour of the solution x(t)
                            ! corresponding to these values of x(0) and r
  ! 
  !==================================================
  ! Modified Euler scheme : a simple Euler step
  ! is used to estimate the derivative as a centered 
  ! (midpoint) difference instead of using the 
  ! leftmost value
  !
  do i=1, N
     t = i*dt
     xt = x + (1. + r*x + x**2) * dt
     x  = x + (2. + r*(x+xt) + x**2 + xt**2)*dt2
     if ( t > tk ) then
         print*, t, x
         tk = tk + ts
     end if 
  end do
  !==================================================
   
  end program






