Hi all, I'm writing a fortran 90 program to simulate the Duffing oscillator. I have it working and have produced a cos wave which I expected for D (the driving force) being set to 0. Here's the 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.0_dp B = -1.0_dp C = 0.0_dp D = 0.0_dp y = 1.0_dp z = 0.0_dp t = 0.0_dp do i=1,999 call rk2(z, y, t, n) n = n + 1.0_dp write(100,*) z, y 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 h = 0.1_dp t = n*h 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 So next I need to set the value for the driving force D to a small value and set it at the unstable equilibrium. What does that mean? How do I set it at the "unstable equilibrium"? Is that just y=0? My second question is about the limit cycle. What does this look like? What should I expect? The question says to adjust the parameters until I observe one, but I don't know what one looks like. I've tried googling it but can't find anything useful. Finally, and most importantly (as it asks me to do this often), it asks me to plot the Poincaré section. To do this, I have to take snap shots of the phase-space diagram at certain time intervals such that t*omega = 2*n*pi, n = 0,1,2,3..... How do I do this, and what should it look like?