\documentclass[12pt]{article}
\newcommand{\dsfrac}[2]{\frac{\displaystyle #1}{\displaystyle #2}}
\author{Mario Pineda-Krch}
\title{Just another tritrophic model}
\begin{document}
\maketitle 
\section{The Hastings \& Powell model}
The non-dimensional version of the Hastings \& Powell (1991) model is given by,
\begin{equation}
\begin{array}{lcl}
\dsfrac{dx}{dt} & = & x(1-x) - \dsfrac{a_1 x}{1+b_1 x}y \\ \\
\dsfrac{dy}{dt} & = & \dsfrac{a_1 x}{1+b_1 x}y - \dsfrac{a_2 y}{1+b_2 y}z - d_1 y \\ \\ 
\dsfrac{dz}{dt} & = & \dsfrac{a_2 y}{1+b_2 y}z - d_2 z 
\end{array}
\end{equation}

Define the nondimensional system in R,
<<>>=
hp = function(time, population, parms){
  x <- population[1]
  y <- population[2]
  z <- population[3]
  with(as.list(parms),{
    dx = x*(1-x) - (a1*x)/(1+b1*x)*y
    dy = (a1*x)/(1+b1*x)*y - (a2*y)/(1+b2*y)*z - d1*y
    dz = (a2*y)/(1+b2*y)*z - d2*z
    out = c(dx, dy, dz)
    list(out)
  })
}
@

I eyeballed the initial conditions from Figure 2 in HP91.
<<>>=
x0 <- c(x=0.75, y=0.15, z=0)
@

Declare the time vector,
<<>>=
time <- seq(0, 5000)
@

Here I am using the same parameters as in Figure 2 in HP91 (which is identical to Figure 2 in KH94),
<<>>=
parms <- c(a1=5.0, b1=2.5, a2=0.1, b2=2.0, d1=0.4, d2=0.01)
@

Solve the system numerically,
<<>>=
require(deSolve)
out <- as.data.frame(ode(x0, time, hp, parms))
@

Look at the structure of the result object,
<<>>=
str(out)
@

and the beginning of the time series,
<<>>=
head(out)
@

\begin{figure}[!h]
\begin{center}
<<fig=TRUE, width=7, height=7>>=
plot(y~x, data=out, cex=.5, pch=19, xlab=expression(italic(x)), ylab=expression(italic(y)))
@
\end{center}
\caption{Phase plane plot of the resource ($x$) and the predator $y$.
}
\end{figure}

\section{Session information}
<<results=tex>>=
toLatex(sessionInfo())
@

\end{document}