Running examples in 'rootSolve-Ex.R' failed
The error most likely occurred in:

> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: plot.steady1D
> ### Title: Plot and Summary Method for steady1D, steady2D and steady3D
> ###   Objects
> ### Aliases: plot.steady1D image.steady2D subset.steady2D image.steady3D
> ###   summary.rootSolve
> ### Keywords: hplot
> 
> ### ** Examples
> 
> ## =======================================================================
> ##  EXAMPLE 1: 1D model, BOD + O2                                
> ## =======================================================================
> ## Biochemical Oxygen Demand (BOD) and oxygen (O2) dynamics
> ## in a river
> 
> #==================#
> # Model equations  #
> #==================#
> O2BOD <- function(t, state, pars) {
+   BOD <- state[1:N]
+   O2  <- state[(N+1):(2*N)]
+ 
+ # BOD dynamics
+   FluxBOD <-  v * c(BOD_0, BOD)  # fluxes due to water transport
+   FluxO2  <-  v * c(O2_0, O2)
+   
+   BODrate <- r*BOD*O2/(O2+10)  # 1-st order consumption, Monod in oxygen
+ 
+ #rate of change = flux gradient - consumption  + reaeration (O2)
+   dBOD         <- -diff(FluxBOD)/dx  - BODrate
+   dO2          <- -diff(FluxO2)/dx   - BODrate + p*(O2sat-O2)
+ 
+   return(list(c(dBOD = dBOD, dO2 = dO2)))
+ }    # END O2BOD
>  
>  
> #==================#
> # Model application#
> #==================#
> # parameters
> dx      <- 100       # grid size, meters
> v       <- 1e2       # velocity, m/day
> x       <- seq(dx/2,10000,by=dx)  # m, distance from river
> N       <- length(x)
> r       <- 0.1       # /day, first-order decay of BOD
> p       <- 0.1       # /day, air-sea exchange rate
> O2sat   <- 300       # mmol/m3 saturated oxygen conc
> O2_0    <- 50        # mmol/m3 riverine oxygen conc
> BOD_0   <- 1500      # mmol/m3 riverine BOD concentration
> 
> # initial guess:
> state <- c(rep(200,N), rep(200,N))
> 
> # running the model
> out   <- steady.1D (y = state, func = O2BOD, parms = NULL,
+                     nspec = 2, pos = TRUE, 
+                     names = c("BOD", "O2"))
> 
> summary(out)
                BOD         O2
Min.       1.444351   2.854498
1st Qu.   14.142193   7.575609
Median   132.765558  98.040933
Mean     377.033511 125.898649
3rd Qu.  696.161549 252.191009
Max.    1430.419020 292.457863
N        100.000000 100.000000
sd       448.960524 115.949342
> 
> # output
> plot(out, grid = x, type = "l", lwd = 2, 
+      ylab = c("mmol/m3", "mmol O2/m3"))
> 
> # observations
> obs <- matrix (ncol = 2, data = c(seq(0, 10000, 2000),
+                                 c(1400, 900,400,100,10,10)))
> 
> colnames(obs) <- c("Distance", "BOD")
> 
> # plot with observations
> plot(out, grid = x, type = "l", lwd = 2, ylab = "mmol/m3", obs = obs, 
+      pch = 16, cex = 1.5)
 ----------- FAILURE REPORT -------------- 
 --- failure: the condition has length > 1 ---
 --- srcref --- 
: 
 --- package (from environment) --- 
rootSolve
 --- call from context --- 
plot.steady1D(out, grid = x, type = "l", lwd = 2, ylab = "mmol/m3", 
    obs = obs, pch = 16, cex = 1.5)
 --- call from argument --- 
if (!class(obs) %in% c("data.frame", "matrix")) stop("'obs' should be either a 'data.frame' or a 'matrix'")
 --- R stacktrace ---
where 1: plot.steady1D(out, grid = x, type = "l", lwd = 2, ylab = "mmol/m3", 
    obs = obs, pch = 16, cex = 1.5)
where 2: plot(out, grid = x, type = "l", lwd = 2, ylab = "mmol/m3", obs = obs, 
    pch = 16, cex = 1.5)

 --- value of length: 2 type: logical ---
[1] FALSE  TRUE
 --- function from context --- 
function (x, ..., which = NULL, grid = NULL, xyswap = FALSE, 
    ask = NULL, obs = NULL, obspar = list(), vertical = FALSE) 
{
    checkX <- function(x) {
        X <- x$y
        if (is.vector(X)) {
            nspec <- attributes(x)$nspec
            if (length(X)%%nspec != 0) 
                stop("length of 'x' should be a multiple of 'nspec' if x is a vector")
            x <- matrix(ncol = nspec, data = X)
        }
        else x <- X
        if (is.null(colnames(x))) 
            colnames(x) <- 1:ncol(x)
        x
    }
    if (!is.null(grid)) 
        if (!is.vector(grid)) 
            stop("'grid' should be a vector")
    preparex <- function(Which, xother, x, xx) {
        ii <- which(xother %in% Which)
        if (length(ii) == length(Which)) 
            xx <- NULL
        if (length(ii) > 0) {
            xnew <- matrix(ncol = length(ii), data = unlist(x[ii + 
                1]))
            colnames(xnew) <- xother[ii]
            xx <- cbind(xx, xnew)
        }
        return(xx)
    }
    xx <- checkX(x)
    nobs <- 0
    if (!is.null(obs)) {
        if (!is.data.frame(obs) & is.list(obs)) {
            Obs <- obs
            obs <- Obs[[1]]
            obs.pos <- matrix(nrow = 1, c(1, nrow(obs)))
            if (!class(obs) %in% c("data.frame", "matrix")) 
                stop("'obs' should be either a 'data.frame' or a 'matrix'")
            if (length(Obs) > 1) 
                for (i in 2:length(Obs)) {
                  obs <- mergeObs(obs, Obs[[i]])
                  obs.pos <- rbind(obs.pos, c(obs.pos[nrow(obs.pos), 
                    2] + 1, nrow(obs)))
                }
            obsname <- colnames(obs)
        }
        else {
            if (is.character(obs[, 1]) | is.factor(obs[, 1])) 
                obs <- convert2wide(obs)
            obsname <- colnames(obs)
            if (!class(obs) %in% c("data.frame", "matrix")) 
                stop("'obs' should be either a 'data.frame' or a 'matrix'")
            obs.pos <- matrix(nrow = 1, c(1, nrow(obs)))
        }
        DD <- duplicated(obsname)
        if (sum(DD) > 0) 
            obs <- mergeObs(obs[, !DD], cbind(obs[, 1], obs[, 
                DD]))
        nobs <- nrow(obs.pos)
    }
    varnames <- c(colnames(xx), names(x)[-1])
    if (is.null(varnames)) 
        varnames <- 1:ncol(xx)
    xother <- names(x)[-1]
    Which <- which
    if (is.null(Which) & is.null(obs)) 
        Which <- 1:ncol(xx)
    else if (is.null(Which)) {
        Which <- which(varnames %in% obsname)
        Which <- varnames[Which]
        if (length(Which) == 0) 
            stop("observed data and model output have no variables in common")
    }
    xx <- preparex(Which, xother, x, xx)
    varnames <- colnames(xx)
    if (length(Which) == 0) 
        stop("nothing to plot")
    xWhich <- selectstvar(Which, varnames)
    np <- length(xWhich)
    ldots <- list(...)
    ndots <- names(ldots)
    ask <- setplotpar(ndots, ldots, np, ask)
    if (ask) {
        oask <- devAskNewPage(TRUE)
        on.exit(devAskNewPage(oask))
    }
    if (!is.null(obs)) {
        ObsWhich <- selectstvar(varnames[xWhich], obsname, NAallowed = TRUE)
        ObsWhich[ObsWhich > ncol(obs)] <- NA
    }
    else ObsWhich <- rep(NA, length(xWhich))
    x2 <- list()
    dots <- list()
    nd <- 0
    nother <- 0
    if (length(ldots) > 0) 
        for (i in 1:length(ldots)) if ("steady1D" %in% class(ldots[[i]])) {
            x2[[nother <- nother + 1]] <- ldots[[i]]
            names(x2)[nother] <- ndots[i]
        }
        else if (is.list(ldots[[i]]) & "steady1D" %in% class(ldots[[i]][[1]])) {
            for (j in 1:length(ldots[[i]])) {
                x2[[nother <- nother + 1]] <- ldots[[i]][[j]]
                names(x2)[nother] <- names(ldots[[i]])[[j]]
            }
        }
        else if (!is.null(ldots[[i]])) {
            dots[[nd <- nd + 1]] <- ldots[[i]]
            names(dots)[nd] <- ndots[i]
        }
    if (nother > 0) {
        for (i in 1:nother) {
            X <- checkX(x2[[i]])
            x2[[i]] <- preparex(Which, xother, x2[[i]], X)
            if (min(dim(x2[[i]]) - dim(xx) == c(0, 0)) == 0) 
                stop(" 'x2' and 'x' are not compatible - dimensions not the same")
            if (min(colnames(x2[[i]]) == varnames) == 0) 
                stop(" 'x2' and 'x' are not compatible - colnames not the same")
        }
    }
    nx <- nother + 1
    if (is.null(grid)) 
        grid <- 1:nrow(xx)
    if (length(grid) != nrow(xx)) 
        stop("length of grid (x-axis) should be = number of rows in 'x$y'")
    plotnames <- c("xlab", "ylab", "xlim", "ylim", "main", "sub", 
        "log", "asp", "ann", "axes", "frame.plot", "panel.first", 
        "panel.last", "cex.lab", "cex.axis", "cex.main")
    ii <- names(dots) %in% plotnames
    dotmain <- dots[ii]
    dotmain <- setdots(dotmain, np)
    dotmain$xlab <- expanddots(dots$xlab, "x", np)
    dotmain$ylab <- expanddots(dots$ylab, "", np)
    dotmain$main <- expanddots(dots$main, varnames[xWhich], np)
    yylim <- expanddotslist(dots$ylim, np)
    xxlim <- expanddotslist(dots$xlim, np)
    ip <- !names(dots) %in% plotnames
    dotpoints <- dots[ip]
    dotpoints <- setdots(dotpoints, nx)
    dotpoints$type <- expanddots(dots$type, "l", nx)
    dotpoints$lty <- expanddots(dots$lty, 1:nx, nx)
    dotpoints$pch <- expanddots(dots$pch, 1:nx, nx)
    dotpoints$col <- expanddots(dots$col, 1:nx, nx)
    dotpoints$bg <- expanddots(dots$bg, 1:nx, nx)
    xyswap <- rep(xyswap, length = np)
    vertical <- rep(vertical, length = np)
    if (nobs > 0) 
        Obspar <- setdots(obspar, nobs)
    for (ip in 1:np) {
        i <- xWhich[ip]
        io <- ObsWhich[ip]
        Dotmain <- extractdots(dotmain, ip)
        Dotpoints <- extractdots(dotpoints, 1)
        Xlog <- Ylog <- FALSE
        if (!is.null(Dotmain$log)) {
            Ylog <- length(grep("y", Dotmain$log))
            Xlog <- length(grep("x", Dotmain$log))
        }
        if (vertical[ip]) {
            xyswap[ip] = TRUE
            Dotmain$axes = FALSE
            Dotmain$xlab = ""
            Dotmain$xaxs = "i"
            Dotmain$yaxs = "i"
        }
        if (!xyswap[ip]) {
            if (is.null(yylim[[ip]])) {
                yrange <- Range(NULL, xx[, i], Ylog)
                if (nother > 0) 
                  for (j in 1:nother) yrange <- Range(yrange, 
                    x2[[j]][, i], Ylog)
                if (!is.na(io)) 
                  yrange <- Range(yrange, obs[, io], Ylog)
                Dotmain$ylim <- yrange
            }
            else Dotmain$ylim <- yylim[[ip]]
            if (Dotmain$xlab == "x" & Dotmain$ylab == "") {
                xl <- Dotmain$ylab
                Dotmain$ylab <- Dotmain$xlab
                Dotmain$xlab <- xl
            }
            if (is.null(xxlim[[ip]])) {
                xrange <- Range(NULL, grid, Xlog)
                if (!is.na(io)) 
                  xrange <- Range(xrange, obs[, 1], Xlog)
                Dotmain$xlim <- xrange
            }
            else Dotmain$xlim <- xxlim[[ip]]
            do.call("plot", c(alist(grid, xx[, i]), Dotmain, 
                Dotpoints))
            if (nother > 0) 
                for (j in 2:nx) do.call("lines", c(alist(grid, 
                  x2[[j - 1]][, i]), extractdots(dotpoints, j)))
            if (!is.na(io)) 
                for (j in 1:nobs) if (length(i.obs <- obs.pos[j, 
                  1]:obs.pos[j, 2]) > 0) 
                  do.call("points", c(alist(obs[i.obs, 1], obs[i.obs, 
                    io]), extractdots(Obspar, j)))
        }
        else {
            if (is.null(xxlim[[ip]])) {
                xrange <- Range(NULL, xx[, i], Xlog)
                if (nother > 0) 
                  for (j in 1:nother) xrange <- Range(xrange, 
                    x2[[j]][, i], Xlog)
                if (!is.na(io)) 
                  xrange <- Range(xrange, obs[, io], Xlog)
                Dotmain$xlim <- xrange
            }
            else Dotmain$xlim <- xxlim[[ip]]
            if (is.null(yylim[[ip]])) {
                yrange <- Range(NULL, range(grid), Ylog)
                if (!is.na(io)) 
                  yrange <- Range(yrange, obs[, 1], Ylog)
                Dotmain$ylim <- rev(yrange)
            }
            else {
                Dotmain$ylim <- yylim[[ip]]
                if (vertical[ip]) 
                  Dotmain$ylim <- Dotmain$ylim[c(2, 1)]
            }
            do.call("plot", c(alist(xx[, i], grid), Dotmain, 
                Dotpoints))
            if (vertical[ip]) {
                abline(h = Dotmain$ylim[2])
                abline(v = Dotmain$xlim[1])
                axis(side = 2)
                axis(side = 3, mgp = c(3, 0.5, 0))
            }
            if (nother > 0) 
                for (j in 2:nx) do.call("lines", c(alist(x2[[j - 
                  1]][, i], grid), extractdots(dotpoints, j)))
            if (!is.na(io)) 
                for (j in 1:nobs) if (length(i.obs <- obs.pos[j, 
                  1]:obs.pos[j, 2]) > 0) 
                  do.call("points", c(alist(obs[i.obs, io], obs[i.obs, 
                    1]), extractdots(Obspar, j)))
        }
    }
}
<bytecode: 0x431dda8>
<environment: namespace:rootSolve>
 --- function search by body ---
Function plot.steady1D in namespace rootSolve has this body.
 ----------- END OF FAILURE REPORT -------------- 
Error in if (!class(obs) %in% c("data.frame", "matrix")) stop("'obs' should be either a 'data.frame' or a 'matrix'") : 
  the condition has length > 1
Calls: plot -> plot.steady1D
Execution halted