library("Matrix") library("optimx") dat <- read.table("min.dat", skip = 3, header = FALSE) nlogL <- function(theta) { k <- exp(theta[1:3]) sigma <- exp(theta[4]) A <- rbind( c(-k[1], k[2]), c( k[1], -(k[2]+k[3])) ) x0 <- c(0,100) sol <- function(t) 100 - sum(expm(A*t)%*%x0) pred <- sapply(dat[,1], sol) -sum(dnorm(dat[,2], mean=pred, sd=sigma, log=TRUE)) } fit <- optimx(c(-2, -2, -2, -2), nlogL, hessian = TRUE, control = list(all.methods = TRUE, save.failures = TRUE, trace = 2)) print(fit)