!
!  This program integrates the one-dimensional flow 
!
!     dx
!     --  = f(x,r)
!     dt
! 
!  with user-supplied functions f(x,r) and f'(x,r)
!  The code also computes the nearest fixed point
!  in the steepest-descent sense and if the FP is 
!  stable also performs a linera stability analysis.
!
!  The code is meant to illustrate the concept of 
!  large-amplitude branches of fixed points, such as 
!  those arising in the ODE implemented here
!
!     f(x,r) = r*x + x**3 - x**5 
!
!  that exhibits a subcritical pitchfork bifurcation 
!  at r = 0 and a saddle-node bifurcation at 
!  r = rc = -1/4 
!
!  Four branches of fixed points (two U and 3 S, xs=0 
!  and two symmetrical large-amplitude branches) exist 
!  in the interval [rx,0]
!
!  The code writes t,x(t),Dx(t) on the file 'ODEint.dat'
!  where Dx(t) is the time evolution of the displacement 
!  between the initial condition and the nearest fixed 
!  point found computed according to a linear-stability
!  analysis.
!
!=========================================================== 
 module global 
   double precision dt,dt2,x,r,t,x0,xf,trel
 end module 
!===========================================================
 program main 
 use global 
 integer NS, is
 double precision TT, f, df
 character (len=1) stab

 TT  = 10.
 dt  = 0.01          ! integration step 
 dt2 = dt/2.d0       ! half step
 is  = 20            ! sampling interval (# steps)
 NS  = int(TT/dt)  
 t   = 0.
 x0  = 0.84
 x   = x0 
 r   = -0.2
!================================================================
! Compute the nearest fixed point (Newton method)

 call fixed_point(x,xf,r)
 trel = 1./df(xf,r)
 stab = 'U'
 if ( trel < 0 ) stab = 'S'
 print*, ' # Nearest fixed point xf (steep. desc.): ', xf
 print*, ' # Linear stability (Unstable/Stable)   : ', stab
 if ( stab == 'S') & 
 print*, ' # Relaxation time                      : ', -trel
 open (1,file='ODEint.dat')
 write(1,*) ' #'
 write(1,*) ' # xf =',xf,' t_rel =', -trel, &
            ' stability: ', stab 
 write(1,*) ' #'
 write(1,*) ' # t           x(t)        (x0-xf)*exp[f''(xf)t]+xf'
 write(1,*) ' #'
 
!================================================================
! Integrate 
!
 do i=1, NS
    t = t + dt
    call one_step_pc
    if ( mod(i,is) == 0 ) call printout
 end do

 end program
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
 subroutine one_step_pc
 use global 
 double precision f,xt
!
! one time step, predictor-corrector method
!
 xt = x + f(x,r) * dt 
 x  = x + (f(x,r) + f(xt,r)) * dt2

end subroutine
!================================================================
 subroutine printout 
 use global 

 write(1,'(3x,3(f9.6,3x))') t, x, (x0-xf)*exp(t/trel)+xf

 end subroutine
!================================================================
 function f(z,rr)
 double precision f,z,rr

 f = rr*z + z**3 - z**5 

 end function
!================================================================
 function df(z,rr)
 double precision df,z,rr

 df = rr + 3*z**2 - 5*z**4 

 end function
!================================================================
 subroutine  fixed_point(xstart,xc,r)
 double precision tol, f, df, xstart, x, xc
 logical tol_not_achieved
 external f, df

 tol = 1.D-10                    !  Specified accuracy
 tol_not_achieved = .true.       !  Switch to end iterations
 x = xstart
 do while ( tol_not_achieved ) 
      xc = x - f(x,r)/df(x,r)
      tol_not_achieved = abs(xc - x) > tol
      x = xc
 end do

 end subroutine

