title: How to plot the Poincaré Section? (Duffing Oscillator)
body: <p>I've written a program that successfully shows a simple limit cycle for the Duffing equation. However, I now need to plot the Poincaré section for this case. </p>
  
  <p>I need to do this by taking snapshots of the Phase-Space diagram at regular time intervals, such that <code>t*omega = 2*pi*n</code>. As I have omega set to 1 for this case, this is just when <code>t = 2*pi*n</code>. I've attempted this, but am not getting the Poincaré section I expect. </p>
  
  <p>Here's my code:</p>
  
  <pre><code>program rungekutta
  implicit none
  integer, parameter :: dp = selected_real_kind(15,300)
  integer :: i, n
  real(kind=dp) z, y, t, A, C, D, B, omega, h
  open(unit=100, file="rungekutta.dat",status='replace')
  n = 0
  !constants
  omega = 1.0_dp
  A = 0.25_dp
  B = 1.0_dp
  C = 0.1_dp
  D = 1.0_dp
  y = 0.0_dp
  z = 0.0_dp
  t = 0.0_dp
  do i=1,1000
  call rk2(z, y, t, n)
  n = n + 1.0_dp
  write(100,*) y, z
  end do
  
  contains
  subroutine rk2(z, y, t, n) !subroutine to implement runge-kutta algorithm
  implicit none
  integer, parameter :: dp = selected_real_kind(15,300)
  integer, intent(inout) :: n
  real(kind=dp) :: k1y, k1z, k2y, k2z, y, z, t, pi
  pi = 4.0*ATAN(1.0)
  h = 0.1_dp
  t = n*2*pi
  k1y = dydt(y,z,t)*h
  k1z = dzdt(y,z,t)*h
  k2z = dzdt(y + (0.5_dp*k1y), z + (0.5_dp*k1z), t + (0.5_dp*h))*h
  k2y = dydt(y, z +(0.5_dp*k1z), t)*h
  y = y + k2y
  z = z + k2z
  end subroutine
  
  !2nd order ODE split into 2 for coupled Runge-Kutta, useful to define 2 
  functions
  function dzdt(y,z,t)
  real(kind=dp) :: y, z, t, dzdt
  dzdt = -A*y**3.0_dp + B*y - C*z + D*sin(omega*t)
  end function
  
  function dydt(y,z,t)
  real(kind=dp) ::  z, dydt, y, t
  dydt = z
  end function
  end program
  </code></pre>
  
  <p>I have also attached an image of what my Poincaré section looks like:</p>
  
  <p><a href="https://i.stack.imgur.com/O3aKB.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/O3aKB.png" alt="Poincaré section"></a>.</p>
  
  <p>This is y on the x axis vs dydt. </p>
  
  <p>And an image of what I'd expect:<a href="https://i.stack.imgur.com/PdERe.gif" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/PdERe.gif" alt="enter image description here"></a></p>