epsilon = .01
theta = fltarr(1001)
theta(0) = !pi / 2.
theta(1) = theta(0)
for n = 1, 999 do begin
  arg = theta(n) - 2. * !pi * floor(theta(n) / (2. * !pi))
  if (arg ge 0.) and (arg lt !pi) then f = 1. else f = -1.
  theta(n + 1) = 2. * theta(n) - f * epsilon^2 * sin(theta(n)) - theta(n - 1)
endfor
plot, theta
end