shipMotn <- function (t, y, pars) {
    ff = approx(pars$ft,pars$fn,t)
    dy = -2*pars$zeta*wn*y[2] - pars$wn*pars$wn*y[1] + ff$y
    dy = dy - pars$zeta_nl*y[1]*y[1]*y[2]
    return(list(c(
      y[2],
      dy
    )))
}

library(deSolve)
wn=28/500*2*pi;
zeta=0.02;
zeta_nl = 0.015;
yini <- c(y1 = -.4, y2 = -.2)
t=(0:6000)/2
ft=t
set.seed(10)
fn=rnorm(length(t))
p=list(wn=wn, zeta=zeta,
    zeta_nl=zeta_nl, ft=ft, fn=fn)
soln <- ode(y = yini, func = shipMotn,
             times = t, parms = p, method='ode23')
plot(soln[1:1000,1],soln[1:1000,2],'l')

library(forecast)
myts <- ts(data=soln[,2],freq=2)
fit <- auto.arima(myts)