Forced oscillation in a pendulum and displacement over time At the old Exploratorium in San Francisco there used to be a Resonant Pendulum.  A weak magnet tied with a string was used to exert a force on a steel plate wrapped around a 200kg concrete cylinder suspended 20m from the ceiling.  We know that for a pendulum without the small angles approximation the equation describing the motion of the mass is:
$$
\ddot{\phi}(t)=-\dfrac{g}{l}\sin\left(\phi(t)\right)
$$
The system begins at rest, and the object is to periodically tug on the string connected to the magnet periodically, using resonance to increase the amplitude of the swing to something noticeable, say a meter on either side of the rest position.
I'm trying to figure out how long it would take to achieve this sort of result for a variety of mass/force relationships.
I assume the force when pulling the string is sinusoidal and the force is zero when the pendulum "pushes" against, instead of being pulled by the string.
Thus, when the amplitude phase is positive, we have
$$
\ddot{\phi}(t)=-\dfrac{g}{l}\sin\left(\phi(t)\right)
$$
and otherwise we have:
$$
\ddot{\phi}(t)=-\dfrac{g}{l}\sin\left(\phi(t)\right)-A\sin(\phi t)
$$
How the heck does one (1) put this into a good differential form, and (2) solve it in something like R (using deSolve)?
Answering my own question (because the comments aren't friendly to code snippets), I think this is the R code I was looking for:
require(deSolve)

force <- 0.5
theta <- 0

pendulum <- function (t, Y, parms) {
    theta <- atan2 (Y['y'], Y['x']) + pi/2
    force <- ifelse(Y['u'] < 0, 0, 0.5)
    with (as.list(Y),
          list(c(u,v,-lam * x + force*cos(theta),-lam * y - 9.8 - force*sin(theta), x^2 + y^2 -1))
          )
}

yini <- c(x = 0, y = -1, u = 0, v = 0, lam = 1)

M <- diag(nrow = 5)
M[5, 5] <- 0

index <- c(2, 2, 1)
times <- seq(from = 0, to = 50, by = 0.01)

out <- radau (y = yini, func = pendulum, parms = NULL, times = times, mass = M, nind = index)

plot(out, type = "l", lwd = 2)
plot(out[, c("x", "y")], type = "l", lwd = 2)

Did I get it right?
 A: OK, the following R code appears to behave well for a wide range of initial and final conditions.  I chose to add damping to the equations because graphing the damped results makes it easy to see whether or what equilibrium state is reached--and whether that is reasonable.  For example, if the magnet force equals the force of gravity, we expect the bob to come to rest at 45 degrees, or at equal values of X and Y.  This is observed.
The R code is as follows:
require(deSolve)

beta <- 0.008                ## Velocity-dependent damping, i.e., wind resistance
gamma <- 0.04                ## Intrinsic damping, i.e., friction of joint
magnet_force <- 0.4          ## Resonant pulling force against bob

pendulum <- function (t, Y, parms) {
    forcing <- ifelse(Y['u'] < 0 || t >= 70, 0, magnet_force)
    with (as.list(Y),
          list(c(u,
                 v,
                 -lam*x - beta*u - gamma*sign(u) + forcing,
                 -lam*y - 9.8 - beta*v - gamma*sign(v),
                 x^2 + y^2 - 1))
          )
}

yini <- c(x = 0, y = -1, u = 0, v = 0, lam = 9.8)

M <- diag(nrow = 5)
M[5, 5] <- 0

index <- c(2, 2, 1)
times <- seq(from = 0, to = 200, by = 0.01)

out <- radau (y = yini, func = pendulum, parms = NULL, times = times, mass = M, nind = index)

plot(out, type = "l", lwd = 2)
plot(out[, c("x", "y")], type = "l", lwd = 2)

A big thanks to Veronica Ciocanel for publishing her paper on pendulums (http://dukespace.lib.duke.edu/dspace/bitstream/handle/10161/5202/dissertationfinal_VeronicaCiocanel.pdf) and including enough other details about the problems she solved that I could tease out what I wanted from her math.
Here is a plot:

