I am trying to fit this very simple 4 species linear Lotka-Volterra competition model to observed data but for some reason when I try the `optim()` function something with regards to `deSolve` seems to fail.

```
# Data
data <- data.frame(Cod = c(0.1966126, 0.1989563, 0.2567677, 0.3158896, 0.4225435, 0.7219856,
                           1.0570824, 0.7266830, 0.6286763, 0.6389475),
                   Herring = c(1.988372, 2.788014, 3.397138, 2.557245, 2.627013, 3.045617, 
                               3.161002, 3.531306, 3.432021, 3.617174),
                   Sprat = c(2.030273, 3.480469, 3.009277, 1.895996, 2.457520, 1.991211, 2.350098,
                             2.118164, 1.693359, 1.869141),
                   Flounder = c(0.4758220, 0.4425532, 0.4185687, 0.4967118, 0.7102515, 0.5733075,
                                0.7404255, 0.5996132, 0.6235977, 0.7187621))
# Model formulation
LLV <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    db1.dt = b1*(r1+a11*b1+a12*b2+a13*b3+a14*b4)
    db2.dt = b2*(r2+a22*b2+a21*b1+a23*b3+a24*b4)
    db3.dt = b3*(r3+a33*b3+a31*b1+a32*b2+a34*b4)
    db4.dt = b4*(r4+a44*b4+a41*b1+a42*b2+a43*b3)
    list(c(db1.dt, db2.dt, db3.dt, db4.dt))
  })
}
# Model input and simulation
# Model input
params <- c(r1 = -0.342085, r2 = 0.6855681, r3 = 2.757769, r4 = 0.9744113,
            a11 = -1.05973762, a12 = 0.09577309, a13 = -0.01915480, a14 = 1.36098939,
            a21 = 0.17533326, a22 = -0.32247342, a23 = 0.03111628, a24 = 0.30212711,
            a31 = 0.5303516, a32 = -0.4869761, a33 = -0.3194882, a34 = -1.5089027,
            a41 = 0.004418133, a42 = 0.163716414, a43 = -0.237873378, a44 = -1.519158802)
ini <- c(b1 = data[1,1], b2 = data[1,2], b3 = data[1,3], b4 = data[1,4])
tmax <- 10
t <- seq(1,tmax,0.1)
# Results and first parameter guess is more or less okay
results <- deSolve::ode(y = ini, times = t, func = LLV, parms = params)
matplot(data, pch = 1)
matplot(x = results[,1], y = results[,-1], type = "l", add = TRUE)

```

Here I proceed and write a function that minimises the residual sum of squares that when included in `optim()` with the above initial parameter guess should produce what I am looking for.

```
min.RSS <- function(data, params) {
  output <- deSolve::ode(y = ini, times = t, func = LLV, parms = params)
  predictions <- exp(output[,-1])
  observations <- data
  return(sum((predictions-observations)^2))
}
result <- optim(par = params, fn = min.RSS, data = data)
fit <- deSolve::ode(y = ini, times = t, func = LLV, parms = result$par)
matplot(x = fit[,1], y = fit[,-1], type = "l", lwd = 3, add = TRUE)

```

Any idea on how to solve this problem will be very much appreciated.