diff --git a/NAMESPACE b/NAMESPACE
index 222bf702..fb719665 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -30,6 +30,7 @@ S3method(getValidNlmixrCtl,fo)
 S3method(getValidNlmixrCtl,foce)
 S3method(getValidNlmixrCtl,focei)
 S3method(getValidNlmixrCtl,foi)
+S3method(getValidNlmixrCtl,nlm)
 S3method(getValidNlmixrCtl,nlme)
 S3method(getValidNlmixrCtl,posthoc)
 S3method(getValidNlmixrCtl,rxSolve)
@@ -57,6 +58,7 @@ S3method(nlmixr2Est,fo)
 S3method(nlmixr2Est,foce)
 S3method(nlmixr2Est,focei)
 S3method(nlmixr2Est,foi)
+S3method(nlmixr2Est,nlm)
 S3method(nlmixr2Est,nlme)
 S3method(nlmixr2Est,output)
 S3method(nlmixr2Est,posthoc)
@@ -130,6 +132,7 @@ S3method(nmObjGetControl,fo)
 S3method(nmObjGetControl,foce)
 S3method(nmObjGetControl,focei)
 S3method(nmObjGetControl,foi)
+S3method(nmObjGetControl,nlm)
 S3method(nmObjGetControl,nlme)
 S3method(nmObjGetControl,posthoc)
 S3method(nmObjGetControl,saem)
@@ -157,6 +160,7 @@ S3method(nmObjGetPredOnly,saem)
 S3method(nmObjGetRxSolve,default)
 S3method(nmObjHandleControlObject,default)
 S3method(nmObjHandleControlObject,foceiControl)
+S3method(nmObjHandleControlObject,nlmControl)
 S3method(nmObjHandleControlObject,nlmeControl)
 S3method(nmObjHandleControlObject,saemControl)
 S3method(nmObjHandleModelObject,default)
@@ -209,9 +213,15 @@ S3method(rxUiGet,foceiThetaS)
 S3method(rxUiGet,getEBEEnv)
 S3method(rxUiGet,getSplitMuModel)
 S3method(rxUiGet,loadPrune)
+S3method(rxUiGet,loadPruneNlm)
 S3method(rxUiGet,loadPruneSaem)
 S3method(rxUiGet,loadPruneSaemPred)
 S3method(rxUiGet,loadPruneSens)
+S3method(rxUiGet,nlmModel0)
+S3method(rxUiGet,nlmParIni)
+S3method(rxUiGet,nlmParName)
+S3method(rxUiGet,nlmParNameFun)
+S3method(rxUiGet,nlmRxModel)
 S3method(rxUiGet,nlmeFD)
 S3method(rxUiGet,nlmeFixedFormula)
 S3method(rxUiGet,nlmeFunction)
@@ -286,6 +296,7 @@ export("%>%")
 export(.foceiPreProcessData)
 export(.nlmixr2estLastPredSimulationInfo)
 export(.nlmixr2objectNameAssign)
+export(.nlmixrNlmFun)
 export(.nlmixrNlmeFun)
 export(.rxGetDVFTransform)
 export(.saemDropMuRefFromModel)
@@ -327,6 +338,7 @@ export(ini)
 export(logit)
 export(lotri)
 export(model)
+export(nlmControl)
 export(nlme)
 export(nlmeControl)
 export(nlmixr)
diff --git a/NEWS.md b/NEWS.md
index ab801ed6..b665c1f1 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,13 @@
 # nlmixr2est (development version)
 
+## New Features
+
+- Algebraic mu referencing has been implemented in `nlme` and `saem`.
+
+- New estimation method "nlm" has been added to estimate population
+  only likelihoods using `stats::nlm` and possibly return a
+  standardized `nlmixr2` fit.
+
 ## Breaking changes
 
 - Removed `fit$saemTransformedData` since it isn't actually used in
@@ -10,9 +18,6 @@
   lag features need to translate the data differently based on
   `rxControl()` options.
 
-## mu referencing
-
-- Algebraic mu referencing has been implemented in `nlme` and `saem`.
 
 ## Bug fixes
 
diff --git a/R/nlm.R b/R/nlm.R
new file mode 100644
index 00000000..dc143102
--- /dev/null
+++ b/R/nlm.R
@@ -0,0 +1,558 @@
+#' nlmixr2 defaults controls for nlm
+#'
+#' @inheritParams stats::nlm
+#' @inheritParams foceiControl
+#' @inheritParams saemControl
+#' @return nlme control object
+#' @details
+#'
+#' Note the covariance is calculated by nlmixr instead of optimHess, so `hessian` is not a possible option
+#'
+#' @export
+#' @author Matthew L. Fidler
+#' @examples
+#'
+#' \donttest{
+#' # A logit regression example with emax model
+#'
+#' dsn <- data.frame(i=1:1000) %>%
+#' dsn$time <- exp(rnorm(1000))
+#' dsn$DV=rbinom(1000,1,exp(-1+dsn$time)/(1+exp(-1+dsn$time)))
+#' dsn$id <- 1 mutate(id=1)
+#'
+#' mod <- function() {
+#'  ini({
+#'    E0 <- 0.5
+#'    Em <- 0.5
+#'    E50 <- 2
+#'    g <- fix(2)
+#'  })
+#'  model({
+#'    v <- E0+Em*time^g/(E50^g+time^g)
+#'    ll(bin) ~ DV * v - log(1 + exp(v))
+#'  })
+#' }
+#'
+#' fit2 <- nlmixr(mod, dsn, est="nlm")
+#'
+#' print(fit2)
+#'
+#' # you can also get the nlm output with fit2$nlm
+#'
+#' fit2$nlm
+#'
+#' # The nlm control has been modified slightly to include
+#' # extra components and name the parameters
+#' }
+nlmControl <- function(typsize = NULL,
+                       fscale = 1, print.level = 2, ndigit = NULL, gradtol = 1e-6,
+                       stepmax = NULL,
+                       steptol = 1e-6, iterlim = 10000, check.analyticals = FALSE,
+                       rxControl=NULL,
+                       optExpression=TRUE, sumProd=FALSE,
+                       returnNlm=FALSE,
+                       addProp = c("combined2", "combined1"),
+                       calcTables=TRUE, compress=TRUE,
+                       covMethod=c("r", "nlm", ""),
+                       adjObf=TRUE, ci=0.95, sigdig=4, sigdigTable=NULL, ...) {
+  checkmate::assertLogical(optExpression, len=1, any.missing=FALSE)
+  checkmate::assertLogical(sumProd, len=1, any.missing=FALSE)
+  checkmate::assertNumeric(stepmax, lower=0, len=1, null.ok=TRUE, any.missing=FALSE)
+  checkmate::assertIntegerish(print.level, lower=1, upper=2, any.missing=FALSE)
+  checkmate::assertNumeric(ndigit, lower=0, len=1, any.missing=FALSE, null.ok=TRUE)
+  checkmate::assertNumeric(gradtol, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertNumeric(steptol, lower=0, len=1, any.missing=FALSE)
+  checkmate::assertIntegerish(iterlim, lower=1, len=1, any.missing=FALSE)
+  checkmate::assertLogical(check.analyticals, len=1, any.missing=FALSE)
+  checkmate::assertLogical(returnNlm, len=1, any.missing=FALSE)
+  checkmate::assertLogical(calcTables, len=1, any.missing=FALSE)
+  checkmate::assertLogical(compress, len=1, any.missing=TRUE)
+  checkmate::assertLogical(adjObf, len=1, any.missing=TRUE)
+  covMethod <- match.arg(covMethod)
+  .xtra <- list(...)
+  .bad <- names(.xtra)
+  .bad <- .bad[!(.bad %in% c("genRxControl"))]
+  if (length(.bad) > 0) {
+    stop("unused argument: ", paste
+    (paste0("'", .bad, "'", sep=""), collapse=", "),
+    call.=FALSE)
+  }
+
+  .genRxControl <- FALSE
+  if (!is.null(.xtra$genRxControl)) {
+    .genRxControl <- .xtra$genRxControl
+  }
+  if (is.null(ndigit)) {
+    ndigit <- sigdig
+  }
+  if (is.null(rxControl)) {
+    if (!is.null(sigdig)) {
+      rxControl <- rxode2::rxControl(sigdig=sigdig)
+    } else {
+      rxControl <- rxode2::rxControl(atol=1e-4, rtol=1e-4)
+    }
+    .genRxControl <- TRUE
+  } else if (inherits(rxControl, "rxControl")) {
+  } else if (is.list(rxControl)) {
+    rxControl <- do.call(rxode2::rxControl, rxControl)
+  } else {
+    stop("solving options 'rxControl' needs to be generated from 'rxode2::rxControl'", call=FALSE)
+  }
+  if (!is.null(sigdig)) {
+    checkmate::assertNumeric(sigdig, lower=1, finite=TRUE, any.missing=TRUE, len=1)
+    if (is.null(sigdigTable)) {
+      sigdigTable <- round(sigdig)
+    }
+  }
+  if (is.null(sigdigTable)) {
+    sigdigTable <- 3
+  }
+  checkmate::assertIntegerish(sigdigTable, lower=1, len=1, any.missing=FALSE)
+
+  .ret <- list(covMethod=covMethod,
+               typsize = typsize,
+               fscale = fscale, print.level = print.level, ndigit=ndigit, gradtol = gradtol,
+               stepmax = stepmax,
+               steptol = steptol, iterlim = iterlim,
+               check.analyticals = check.analyticals,
+               optExpression=optExpression,
+               sumProd=sumProd,
+               rxControl=rxControl,
+               returnNlm=returnNlm, addProp=addProp, calcTables=calcTables,
+               compress=compress,
+               ci=ci, sigdig=sigdig, sigdigTable=sigdigTable,
+               genRxControl=.genRxControl)
+  class(.ret) <- "nlmControl"
+  .ret
+}
+
+#' Get the nlm family control
+#'
+#' @param env nlme optimization environment
+#' @param ... Other arguments
+#' @return Nothing, called for side effects
+#' @author Matthew L. Fidler
+#' @noRd
+.nlmFamilyControl <- function(env, ...) {
+  .ui <- env$ui
+  .control <- env$control
+  if (is.null(.control)) {
+    .control <- nlmixr2est::nlmControl()
+  }
+  if (!inherits(.control, "nlmControl")){
+    .control <- do.call(nlmixr2est::nlmControl, .control)
+  }
+  assign("control", .control, envir=.ui)
+}
+
+
+#' @rdname nmObjHandleControlObject
+#' @export
+nmObjHandleControlObject.nlmControl <- function(control, env) {
+  assign("nlmControl", control, envir=env)
+}
+
+#' @rdname nmObjGetControl
+#' @export
+nmObjGetControl.nlm <- function(x, ...) {
+  .env <- x[[1]]
+  if (exists("nlmControl", .env)) {
+    .control <- get("nlmControl", .env)
+    if (inherits(.control, "nlmControl")) return(.control)
+  }
+  if (exists("control", .env)) {
+    .control <- get("control", .env)
+    if (inherits(.control, "nlmControl")) return(.control)
+  }
+  stop("cannot find nlm related control object", call.=FALSE)
+}
+
+#' @rdname getValidNlmixrControl
+#' @export
+getValidNlmixrCtl.nlm <- function(control) {
+  .ctl <- control[[1]]
+  if (is.null(.ctl)) .ctl <- nlmControl()
+  if (is.null(attr(.ctl, "class")) && is(.ctl, "list")) .ctl <- do.call("nlmControl", .ctl)
+  if (!inherits(.ctl, "nlmControl")) {
+    .minfo("invalid control for `est=\"nlm\"`, using default")
+    .ctl <- nlmControl()
+  } else {
+    .ctl <- do.call(nlmControl, .ctl)
+  }
+  .ctl
+}
+
+.nlmEnv <- new.env(parent=emptyenv())
+
+#' A surrogate function for nlm to call for ode solving
+#'
+#' @param dv The observations for the `nlm` function
+#' @param pars Parameters that will be estimated
+#' @return Predictions
+#' @details
+#' This is an internal function and should not be called directly.
+#' @author Matthew L. Fidler
+#' @keywords internal
+#' @export
+.nlmixrNlmFun <- function(pars) {
+  .retF <- do.call(rxode2::rxSolve,
+                   c(list(object=.nlmEnv$model,
+                          params=.nlmEnv$parTrans(pars),
+                          events=.nlmEnv$data),
+                     .nlmEnv$rxControl))
+  -sum(.retF$rx_pred_)
+}
+
+
+#'@export
+rxUiGet.nlmModel0 <- function(x, ...) {
+  .ui <- rxode2::rxUiDecompress(x[[1]])
+  .predDf <- .ui$predDf
+  .save <- .predDf
+  .predDf[.predDf$distribution == "norm", "distribution"] <- "dnorm"
+  assign("predDf", .predDf, envir=.ui)
+  on.exit(assign("predDf", .save, envir=.ui))
+  .ui$foceiModel0ll
+}
+
+#' Load the saem model into symengine
+#'
+#' @param x rxode2 UI object
+#' @return String for loading into symengine
+#' @author Matthew L. Fidler
+#' @noRd
+.nlmPrune <- function(x) {
+  .x <- x[[1]]
+  .x <- .x$nlmModel0[[-1]]
+  .env <- new.env(parent = emptyenv())
+  .env$.if <- NULL
+  .env$.def1 <- NULL
+  .malert("pruning branches ({.code if}/{.code else}) of nlm model...")
+  .ret <- rxode2::.rxPrune(.x, envir = .env)
+  .mv <- rxode2::rxModelVars(.ret)
+  ## Need to convert to a function
+  if (rxode2::.rxIsLinCmt() == 1L) {
+    .vars <- c(.mv$params, .mv$lhs, .mv$slhs)
+    .mv <- rxode2::.rxLinCmtGen(length(.mv$state), .vars)
+  }
+  .msuccess("done")
+  rxode2::rxNorm(.mv)
+}
+
+#' @export
+rxUiGet.loadPruneNlm <- function(x, ...) {
+  .loadSymengine(.nlmPrune(x), promoteLinSens = FALSE)
+}
+
+#' @export
+rxUiGet.nlmRxModel <- function(x, ...) {
+  .s <- rxUiGet.loadPruneNlm(x, ...)
+  .prd <- get("rx_pred_", envir = .s)
+  .prd <- paste0("rx_pred_=", rxode2::rxFromSE(.prd))
+  ## .lhs0 <- .s$..lhs0
+  ## if (is.null(.lhs0)) .lhs0 <- ""
+  .ddt <- .s$..ddt
+  if (is.null(.ddt)) .ddt <- ""
+  .ret <- paste(c(
+    #.s$..stateInfo["state"],
+    #.lhs0,
+    .ddt,
+    .prd,
+    #.s$..stateInfo["statef"],
+    #.s$..stateInfo["dvid"],
+    ""
+  ), collapse = "\n")
+  .sumProd <- rxode2::rxGetControl(x[[1]], "sumProd", FALSE)
+  .optExpression <- rxode2::rxGetControl(x[[1]], "optExpression", TRUE)
+  if (.sumProd) {
+    .malert("stabilizing round off errors in nlm model...")
+    .ret <- rxode2::rxSumProdModel(.ret)
+    .msuccess("done")
+  }
+  if (.optExpression) {
+    .ret <- rxode2::rxOptExpr(.ret, "nlm model")
+    .msuccess("done")
+  }
+  paste(c(rxUiGet.foceiParams(x, ...), rxUiGet.foceiCmtPreModel(x, ...),
+          .ret, .foceiToCmtLinesAndDvid(x[[1]])), collapse="\n")
+}
+
+#' @export
+rxUiGet.nlmParNameFun <- function(x, ...) {
+  .ui <- x[[1]]
+  .iniDf <- .ui$iniDf
+  .env <- new.env(parent=emptyenv())
+  .env$i <- 1
+  eval(str2lang(
+    paste0("function(p) {c(",
+           paste(vapply(seq_along(.iniDf$ntheta), function(t) {
+             if (.iniDf$fix[t]) {
+               paste0("'THETA[", t, "]'=", .iniDf$est[t])
+             } else {
+               if (!is.na(.iniDf$err[t]) &&
+                     !is.finite(.iniDf$upper[t]) &&
+                      .iniDf$lower[t] == 0) {
+                 .ret <- paste0("'THETA[", t, "]'=abs(p[", .env$i, "])")
+               } else {
+                 .ret <- paste0("'THETA[", t, "]'=p[", .env$i, "]")
+                 .env$i <- .env$i + 1
+               }
+               .ret
+             }
+           }, character(1), USE.NAMES=FALSE), collapse=","), ")}")))
+}
+
+#' @export
+rxUiGet.nlmParIni <- function(x, ...) {
+  .ui <- x[[1]]
+  .ui$iniDf$est[!.ui$iniDf$fix]
+}
+
+#' @export
+rxUiGet.nlmParName <- function(x, ...) {
+  .ui <- x[[1]]
+  .ui$iniDf$name[!.ui$iniDf$fix]
+}
+
+#' Setup the data for nlm estimation
+#'
+#' @param dataSav Formatted Data
+#' @return Nothing, called for side effects
+#' @author Matthew L. Fidler
+#' @noRd
+.nlmFitDataSetup <- function(dataSav) {
+  .dsAll <- dataSav[dataSav$EVID != 2, ] # Drop EVID=2 for estimation
+  if (any(names(.dsAll) == "CENS")) {
+    if (!all(.dsAll$CENS == 0)) {
+      stop("'nlm' does not work with censored data", call. =FALSE)
+    }
+  }
+  .nlmEnv$data <- rxode2::etTrans(.dsAll, .nlmEnv$model)
+}
+
+.nlmFitModel <- function(ui, dataSav) {
+  .nlmEnv$model <- rxode2::rxode2(ui$nlmRxModel)
+  .nlmFitDataSetup(dataSav)
+  .nlmEnv$rxControl <- rxode2::rxGetControl(ui, "rxControl", rxode2::rxControl())
+  .nlmEnv$rxControl$returnType <- 2L # use data.frame output
+  .nlmEnv$parTrans <- ui$nlmParNameFun
+  .ctl <- ui$control
+  class(.ctl) <- NULL
+  .p <- ui$nlmParIni
+  .typsize <- .ctl$typsize
+  if (is.null(.typsize)) {
+    .typsize <- rep(1, length(.p))
+  } else if (length(.typsize) == 1L) {
+    .typsize <- rep(.typsize, length(.p))
+  } else {
+    stop("'typsize' needs to match the number of estimated parameters (or equal 1)", call.=FALSE)
+  }
+  .stepmax <- .ctl$stepmax
+  if (is.null(.stepmax)) {
+    .stepmax <- max(1000 * sqrt(sum((.p/.typsize)^2)), 1000)
+  }
+  .ret <- eval(bquote(stats::nlm(
+    f=.(nlmixr2est::.nlmixrNlmFun),
+    p=.(.p),
+    hessian=.(.ctl$covMethod == "nlm"),
+    typsize=.(.typsize),
+    fscale=.(.ctl$fscale),
+    print.level=.(.ctl$print.level),
+    ndigit=.(.ctl$ndigit),
+    gradtol=.(.ctl$gradtol),
+    stepmax=.(.stepmax),
+    steptol = .(.ctl$steptol),
+    iterlim = .(.ctl$iterlim),
+    check.analyticals = .(.ctl$check.analyticals)
+  )))
+  # be nice and name items
+  .name <- ui$nlmParName
+  names(.ret$estimate) <- .name
+  names(.ret$gradient) <- .name
+  if (any(.ctl$covMethod == c("r", "nlm"))) {
+    .malert("calculating covariance")
+    if (.ctl$covMethod != "nlm") {
+      .p <- setNames(.ret$estimate, NULL)
+      .hess <- nlmixr2Hess(.p, nlmixr2est::.nlmixrNlmFun)
+      .ret$hessian <- .hess
+    }
+    dimnames(.ret$hessian) <- list(.name, .name)
+    .hess <- .ret$hessian
+
+    # r matrix
+    .r <- 0.5 * .hess
+    .ch <- try(cholSE(.r), silent = TRUE)
+    .covType <- "r"
+    if (inherits(.ch, "try-error")) {
+      .r2 <- .r %*% .r
+      .r2 <- try(sqrtm(.r2), silent=TRUE)
+      .covType <- "|r|"
+      if (!inherits(.r2, "try-error")) {
+        .ch <- try(cholSE(.r), silent=TRUE)
+        if (inherits(.ch, "try-error")) {
+          .r2 <- .ch # switch to nearPD
+        }
+      }
+      if (inherits(.r2, "try-error")) {
+        .covType <- "r+"
+        .r2 <- try(nmNearPD(.r), silent=TRUE)
+        if (!inherits(.r2, "try-error")) {
+          .ch <- try(cholSE(.r), silent=TRUE)
+        }
+      } else {
+        .ch <- try(cholSE(.r), silent=TRUE)
+      }
+    }
+    if (!inherits(.ch, "try-error")) {
+      .rinv <- rxode2::rxInv(.ch)
+      .rinv <- .rinv %*% t(.rinv)
+      .cov <- 2*.rinv
+      dimnames(.cov) <- list(.name, .name)
+      dimnames(.rinv) <- list(.name, .name)
+      .ret$covMethod <- .covType
+      if (.ctl$covMethod == "nlm") {
+        .ret$covMethod <- paste0(.covType, " (nlm)")
+      } else {
+        .ret$covMethod <- .covType
+      }
+      .ret$cov <- .cov
+
+    } else {
+      .ret$covMethod <- "failed"
+    }
+    dimnames(.r) <- list(.name, .name)
+    .ret$r <- .r
+    .msuccess("done")
+  }
+  .ret
+}
+#' Get the full theta for nlm methods
+#'
+#' @param nlm enhanced nlm return
+#' @param ui ui object
+#' @return named theta matrix
+#' @author Matthew L. Fidler
+#' @noRd
+.nlmGetTheta <- function(nlm, ui) {
+  .iniDf <- ui$iniDf
+  setNames(vapply(seq_along(.iniDf$name),
+         function(i) {
+           if (.iniDf$fix[i]) {
+             return(.iniDf$est[i])
+           } else {
+             return(nlm$estimate[.iniDf$name[i]])
+           }
+         }, double(1), USE.NAMES=FALSE),
+         .iniDf$name)
+}
+
+.nlmControlToFoceiControl <- function(env, assign=TRUE) {
+  .nlmControl <- env$nlmControl
+  .ui <- env$ui
+  .foceiControl <- foceiControl(rxControl=env$nlmControl$rxControl,
+                                maxOuterIterations=0L,
+                                maxInnerIterations=0L,
+                                covMethod=0L,
+                                sumProd=.nlmControl$sumProd,
+                                optExpression=.nlmControl$optExpression,
+                                scaleTo=0,
+                                calcTables=.nlmControl$calcTables,
+                                addProp=.nlmControl$addProp,
+                                #skipCov=.ui$foceiSkipCov,
+                                interaction=0L,
+                                compress=.nlmControl$compress,
+                                ci=.nlmControl$ci,
+                                sigdigTable=.nlmControl$sigdigTable)
+  if (assign) env$control <- .foceiControl
+  .foceiControl
+}
+
+
+.nlmFamilyFit <- function(env, ...) {
+  .ui <- env$ui
+  .control <- .ui$control
+  .data <- env$data
+  .ret <- new.env(parent=emptyenv())
+  # The environment needs:
+  # - table for table options
+  # - $origData -- Original Data
+  # - $dataSav -- Processed data from .foceiPreProcessData
+  # - $idLvl -- Level information for ID factor added
+  # - $covLvl -- Level information for items to convert to factor
+  # - $ui for ui fullTheta Full theta information
+  # - $etaObf data frame with ID, etas and OBJI
+  # - $cov For covariance
+  # - $covMethod for the method of calculating the covariance
+  # - $adjObf Should the objective function value be adjusted
+  # - $objective objective function value
+  # - $extra Extra print information
+  # - $method Estimation method (for printing)
+  # - $omega Omega matrix
+  # - $theta Is a theta data frame
+  # - $model a list of model information for table generation.  Needs a `predOnly` model
+  # - $message Message for display
+  # - $est estimation method
+  # - $ofvType (optional) tells the type of ofv is currently being used
+  # When running the focei problem to create the nlmixr object, you also need a
+  #  foceiControl object
+  .ret$table <- env$table
+  .foceiPreProcessData(.data, .ret, .ui, .control$rxControl)
+  .nlm <- .collectWarn(.nlmFitModel(.ui, .ret$dataSav), lst = TRUE)
+  .ret$nlm <- .nlm[[1]]
+  .ret$message <- NULL
+  if (rxode2::rxGetControl(.ui, "returnNlm", FALSE)) {
+    return(.ret$nlm)
+  }
+  if (.ret$nlm$code == 1) {
+    .ret$message <- "relative gradient is close to zero, current iterate is probably solution"
+  } else if (.ret$nlm$code == 2) {
+    .ret$message <- "successive iterates within tolerance, current iterate is probably solution"
+  } else if (.ret$nlm$code == 3) {
+    .ret$message <- c("last global step failed to locate a point lower than 'estimate'",
+                      "either 'estimate' is an approximate local minimum of the function or 'steptol' is too small")
+  } else if (.ret$nlm$code == 4) {
+    .ret$message <- "iteration limit exceeded"
+  } else if (.ret$nlm$code == 5) {
+    .ret$message <- c("maximum step size 'stepmax' exceeded five consecutive times",
+                      "either the function is unbounded below, becomes asymptotic to a finite value from above in some direction or 'stepmax' is too small")
+  } else {
+    .ret$message <- ""
+  }
+  .ret$ui <- .ui
+  .ret$adjObf <- rxode2::rxGetControl(.ui, "adjObf", TRUE)
+  .ret$fullTheta <- .nlmGetTheta(.ret$nlm, .ui)
+  .ret$cov <- .ret$nlm$cov
+  .ret$covMethod <- .ret$nlm$covMethod
+  #.ret$etaMat <- NULL
+  #.ret$etaObf <- NULL
+  #.ret$omega <- NULL
+  .ret$control <- .control
+  .ret$extra <- ""
+  .nlmixr2FitUpdateParams(.ret)
+  nmObjHandleControlObject(.ret$control, .ret)
+  if (exists("control", .ui)) {
+    rm(list="control", envir=.ui)
+  }
+  .ret$est <- "nlm"
+  # There is no parameter history for nlme
+  .ret$objective <- -2 * as.numeric(.ret$nlm$minimum)
+  .ret$model <- .ui$ebe
+  .ret$ofvType <- "nlm"
+  .nlmControlToFoceiControl(.ret)
+  .ret$theta <- .ret$ui$saemThetaDataFrame
+  .ret <- nlmixr2CreateOutputFromUi(.ret$ui, data=.ret$origData, control=.ret$control, table=.ret$table, env=.ret, est="nlm")
+  .env <- .ret$env
+  .env$method <- "nlm"
+  .ret
+}
+
+#' @rdname nlmixr2Est
+#' @export
+nlmixr2Est.nlm <- function(env, ...) {
+  .ui <- env$ui
+  rxode2::assertRxUiPopulationOnly(.ui, " for the estimation routine 'nlm', try 'focei'", .var.name=.ui$modelName)
+  rxode2::assertRxUiRandomOnIdOnly(.ui, " for the estimation routine 'nlm'", .var.name=.ui$modelName)
+  .nlmFamilyControl(env, ...)
+  on.exit({if (exists("control", envir=.ui)) rm("control", envir=.ui)}, add=TRUE)
+  .nlmFamilyFit(env,  ...)
+}
diff --git a/man/dot-nlmixrNlmFun.Rd b/man/dot-nlmixrNlmFun.Rd
new file mode 100644
index 00000000..f8a7f61b
--- /dev/null
+++ b/man/dot-nlmixrNlmFun.Rd
@@ -0,0 +1,26 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/nlm.R
+\name{.nlmixrNlmFun}
+\alias{.nlmixrNlmFun}
+\title{A surrogate function for nlm to call for ode solving}
+\usage{
+.nlmixrNlmFun(pars)
+}
+\arguments{
+\item{pars}{Parameters that will be estimated}
+
+\item{dv}{The observations for the `nlm` function}
+}
+\value{
+Predictions
+}
+\description{
+A surrogate function for nlm to call for ode solving
+}
+\details{
+This is an internal function and should not be called directly.
+}
+\author{
+Matthew L. Fidler
+}
+\keyword{internal}
diff --git a/man/getValidNlmixrControl.Rd b/man/getValidNlmixrControl.Rd
index 6ab8c611..88c90e47 100644
--- a/man/getValidNlmixrControl.Rd
+++ b/man/getValidNlmixrControl.Rd
@@ -1,6 +1,7 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/sharedControl.R
-\name{getValidNlmixrControl}
+% Please edit documentation in R/nlm.R, R/sharedControl.R
+\name{getValidNlmixrCtl.nlm}
+\alias{getValidNlmixrCtl.nlm}
 \alias{getValidNlmixrControl}
 \alias{getValidNlmixrCtl}
 \alias{getValidNlmixrCtl.focei}
@@ -17,6 +18,8 @@
 \alias{getValidNlmixrCtl.default}
 \title{Get valid nlmixr control object}
 \usage{
+\method{getValidNlmixrCtl}{nlm}(control)
+
 getValidNlmixrControl(control, est)
 
 getValidNlmixrCtl(control)
diff --git a/man/nlmControl.Rd b/man/nlmControl.Rd
new file mode 100644
index 00000000..e879b7fa
--- /dev/null
+++ b/man/nlmControl.Rd
@@ -0,0 +1,205 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/nlm.R
+\name{nlmControl}
+\alias{nlmControl}
+\title{nlmixr2 defaults controls for nlm}
+\usage{
+nlmControl(
+  typsize = NULL,
+  fscale = 1,
+  print.level = 2,
+  ndigit = NULL,
+  gradtol = 1e-06,
+  stepmax = NULL,
+  steptol = 1e-06,
+  iterlim = 10000,
+  check.analyticals = FALSE,
+  rxControl = NULL,
+  optExpression = TRUE,
+  sumProd = FALSE,
+  returnNlm = FALSE,
+  addProp = c("combined2", "combined1"),
+  calcTables = TRUE,
+  compress = TRUE,
+  covMethod = c("r", "nlm", ""),
+  adjObf = TRUE,
+  ci = 0.95,
+  sigdig = 4,
+  sigdigTable = NULL,
+  ...
+)
+}
+\arguments{
+\item{typsize}{an estimate of the size of each parameter
+    at the minimum.}
+
+\item{fscale}{an estimate of the size of \code{f} at the minimum.}
+
+\item{print.level}{this argument determines the level of printing
+    which is done during the minimization process.  The default
+    value of \code{0} means that no printing occurs, a value of \code{1}
+    means that initial and final details are printed and a value
+    of 2 means that full tracing information is printed.}
+
+\item{ndigit}{the number of significant digits in the function \code{f}.}
+
+\item{gradtol}{a positive scalar giving the tolerance at which the
+    scaled gradient is considered close enough to zero to
+    terminate the algorithm.  The scaled gradient is a
+    measure of the relative change in \code{f} in each direction
+    \code{p[i]} divided by the relative change in \code{p[i]}.}
+
+\item{stepmax}{a positive scalar which gives the maximum allowable
+    scaled step length.  \code{stepmax} is used to prevent steps which
+    would cause the optimization function to overflow, to prevent the
+    algorithm from leaving the area of interest in parameter space, or to
+    detect divergence in the algorithm. \code{stepmax} would be chosen
+    small enough to prevent the first two of these occurrences, but should
+    be larger than any anticipated reasonable step.}
+
+\item{steptol}{A positive scalar providing the minimum allowable
+    relative step length.}
+
+\item{iterlim}{a positive integer specifying the maximum number of
+    iterations to be performed before the program is terminated.}
+
+\item{check.analyticals}{a logical scalar specifying whether the
+    analytic gradients and Hessians, if they are supplied, should be
+    checked against numerical derivatives at the initial parameter
+    values. This can help detect incorrectly formulated gradients or
+    Hessians.}
+
+\item{rxControl}{`rxode2` ODE solving options during fitting, created with `rxControl()`}
+
+\item{optExpression}{Optimize the rxode2 expression to speed up
+calculation. By default this is turned on.}
+
+\item{sumProd}{Is a boolean indicating if the model should change
+multiplication to high precision multiplication and sums to
+high precision sums using the PreciseSums package.  By default
+this is \code{FALSE}.}
+
+\item{addProp}{specifies the type of additive plus proportional
+  errors, the one where standard deviations add (combined1) or the
+  type where the variances add (combined2).
+
+The combined1 error type can be described by the following equation:
+
+  y = f + (a + b*f^c)*err
+
+The combined2 error model can be described by the following equation:
+
+ y = f + sqrt(a^2 + b^2*(f^c)^2)*err
+
+ Where:
+
+ - y represents the observed value
+
+ - f represents the predicted value
+
+ - a  is the additive standard deviation
+
+ - b is the proportional/power standard deviation
+
+ - c is the power exponent (in the proportional case c=1)}
+
+\item{calcTables}{This boolean is to determine if the foceiFit
+will calculate tables. By default this is \code{TRUE}}
+
+\item{compress}{Should the object have compressed items}
+
+\item{covMethod}{Method for calculating covariance.  In this
+    discussion, R is the Hessian matrix of the objective
+    function. The S matrix is the sum of individual
+    gradient cross-product (evaluated at the individual empirical
+    Bayes estimates).
+
+\itemize{
+
+ \item "\code{r,s}" Uses the sandwich matrix to calculate the
+ covariance, that is: \code{solve(R) \%*\% S \%*\% solve(R)}
+
+ \item "\code{r}" Uses the Hessian matrix to calculate the
+ covariance as \code{2 \%*\% solve(R)}
+
+ \item "\code{s}" Uses the cross-product matrix to calculate the
+ covariance as \code{4 \%*\% solve(S)}
+
+ \item "" Does not calculate the covariance step.
+}}
+
+\item{adjObf}{is a boolean to indicate if the objective function
+should be adjusted to be closer to NONMEM's default objective
+function.  By default this is \code{TRUE}}
+
+\item{ci}{Confidence level for some tables.  By default this is
+0.95 or 95\% confidence.}
+
+\item{sigdig}{Optimization significant digits. This controls:
+
+\itemize{
+
+ \item The tolerance of the inner and outer optimization is \code{10^-sigdig}
+
+ \item The tolerance of the ODE solvers is
+ \code{0.5*10^(-sigdig-2)}; For the sensitivity equations and
+ steady-state solutions the default is \code{0.5*10^(-sigdig-1.5)}
+ (sensitivity changes only applicable for liblsoda)
+
+ \item The tolerance of the boundary check is \code{5 * 10 ^ (-sigdig + 1)}
+
+}}
+
+\item{sigdigTable}{Significant digits in the final output table.
+If not specified, then it matches the significant digits in the
+`sigdig` optimization algorithm.  If `sigdig` is NULL, use 3.}
+
+\item{...}{additional arguments to be passed to \code{f}.}
+}
+\value{
+nlme control object
+}
+\description{
+nlmixr2 defaults controls for nlm
+}
+\details{
+Note the covariance is calculated by nlmixr instead of optimHess, so `hessian` is not a possible option
+}
+\examples{
+
+\donttest{
+# A logit regression example with emax model
+
+dsn <- data.frame(i=1:1000) \%>\%
+dsn$time <- exp(rnorm(1000))
+dsn$DV=rbinom(1000,1,exp(-1+dsn$time)/(1+exp(-1+dsn$time)))
+dsn$id <- 1 mutate(id=1)
+
+mod <- function() {
+ ini({
+   E0 <- 0.5
+   Em <- 0.5
+   E50 <- 2
+   g <- fix(2)
+ })
+ model({
+   v <- E0+Em*time^g/(E50^g+time^g)
+   ll(bin) ~ DV * v - log(1 + exp(v))
+ })
+}
+
+fit2 <- nlmixr(mod, dsn, est="nlm")
+
+print(fit2)
+
+# you can also get the nlm output with fit2$nlm
+
+fit2$nlm
+
+# The nlm control has been modified slightly to include
+# extra components and name the parameters
+}
+}
+\author{
+Matthew L. Fidler
+}
diff --git a/man/nlmixr2Est.Rd b/man/nlmixr2Est.Rd
index 0f9138e9..35fa23d2 100644
--- a/man/nlmixr2Est.Rd
+++ b/man/nlmixr2Est.Rd
@@ -1,6 +1,6 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/focei.R, R/nlme.R, R/nlmixr2Est.R, R/rxsolve.R,
-%   R/saem.R
+% Please edit documentation in R/focei.R, R/nlm.R, R/nlme.R, R/nlmixr2Est.R,
+%   R/rxsolve.R, R/saem.R
 \name{nlmixr2Est.focei}
 \alias{nlmixr2Est.focei}
 \alias{nlmixr2Est.foce}
@@ -8,6 +8,7 @@
 \alias{nlmixr2Est.foi}
 \alias{nlmixr2Est.fo}
 \alias{nlmixr2Est.output}
+\alias{nlmixr2Est.nlm}
 \alias{nlmixr2Est.nlme}
 \alias{nlmixr2Est}
 \alias{nlmixr2Est.default}
@@ -30,6 +31,8 @@
 
 \method{nlmixr2Est}{output}(env, ...)
 
+\method{nlmixr2Est}{nlm}(env, ...)
+
 \method{nlmixr2Est}{nlme}(env, ...)
 
 nlmixr2Est(env, ...)
diff --git a/man/nmObjGetControl.Rd b/man/nmObjGetControl.Rd
index db6428cb..6264cb76 100644
--- a/man/nmObjGetControl.Rd
+++ b/man/nmObjGetControl.Rd
@@ -1,6 +1,7 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/nlme.R, R/nmObjHandle.R
-\name{nmObjGetControl.nlme}
+% Please edit documentation in R/nlm.R, R/nlme.R, R/nmObjHandle.R
+\name{nmObjGetControl.nlm}
+\alias{nmObjGetControl.nlm}
 \alias{nmObjGetControl.nlme}
 \alias{nmObjGetControl}
 \alias{nmObjGetControl.focei}
@@ -12,6 +13,8 @@
 \alias{nmObjGetControl.default}
 \title{Get control object from fit}
 \usage{
+\method{nmObjGetControl}{nlm}(x, ...)
+
 \method{nmObjGetControl}{nlme}(x, ...)
 
 nmObjGetControl(x, ...)
diff --git a/man/nmObjHandleControlObject.Rd b/man/nmObjHandleControlObject.Rd
index f9e4dc41..9941d87a 100644
--- a/man/nmObjHandleControlObject.Rd
+++ b/man/nmObjHandleControlObject.Rd
@@ -1,6 +1,7 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/nlme.R, R/nmObjHandle.R
-\name{nmObjHandleControlObject.nlmeControl}
+% Please edit documentation in R/nlm.R, R/nlme.R, R/nmObjHandle.R
+\name{nmObjHandleControlObject.nlmControl}
+\alias{nmObjHandleControlObject.nlmControl}
 \alias{nmObjHandleControlObject.nlmeControl}
 \alias{nmObjHandleControlObject}
 \alias{nmObjHandleControlObject.foceiControl}
@@ -8,6 +9,8 @@
 \alias{nmObjHandleControlObject.default}
 \title{Handle the control object}
 \usage{
+\method{nmObjHandleControlObject}{nlmControl}(control, env)
+
 \method{nmObjHandleControlObject}{nlmeControl}(control, env)
 
 nmObjHandleControlObject(control, env)
diff --git a/man/reexports.Rd b/man/reexports.Rd
index 2b382cb2..5560e195 100644
--- a/man/reexports.Rd
+++ b/man/reexports.Rd
@@ -5,6 +5,7 @@
 \alias{reexports}
 \alias{\%>\%}
 \alias{rxode2}
+\alias{as.rxUi}
 \alias{rxode}
 \alias{RxODE}
 \alias{ini}
@@ -75,6 +76,6 @@ below to see their documentation.
 
   \item{nlme}{\code{\link[nlme]{ACF}}, \code{\link[nlme]{augPred}}, \code{\link[nlme]{fixed.effects}}, \code{\link[nlme:fixed.effects]{fixef}}, \code{\link[nlme]{getData}}, \code{\link[nlme]{getVarCov}}, \code{\link[nlme]{groupedData}}, \code{\link[nlme]{nlme}}, \code{\link[nlme]{pdBlocked}}, \code{\link[nlme]{pdCompSymm}}, \code{\link[nlme]{pdConstruct}}, \code{\link[nlme]{pdDiag}}, \code{\link[nlme]{pdFactor}}, \code{\link[nlme]{pdIdent}}, \code{\link[nlme]{pdLogChol}}, \code{\link[nlme]{pdMat}}, \code{\link[nlme]{pdMatrix}}, \code{\link[nlme]{pdNatural}}, \code{\link[nlme]{pdSymm}}, \code{\link[nlme]{random.effects}}, \code{\link[nlme:random.effects]{ranef}}, \code{\link[nlme]{reStruct}}, \code{\link[nlme]{varComb}}, \code{\link[nlme]{varConstPower}}, \code{\link[nlme]{VarCorr}}, \code{\link[nlme]{varExp}}, \code{\link[nlme]{varFixed}}, \code{\link[nlme]{varFunc}}, \code{\link[nlme]{varIdent}}, \code{\link[nlme]{varPower}}, \code{\link[nlme]{varWeights}}}
 
-  \item{rxode2}{\code{\link[rxode2:reexports]{add.dosing}}, \code{\link[rxode2:reexports]{add.sampling}}, \code{\link[rxode2:reexports]{et}}, \code{\link[rxode2:reexports]{et}}, \code{\link[rxode2:reexports]{eventTable}}, \code{\link[rxode2:logit]{expit}}, \code{\link[rxode2:stat_amt]{geom_amt}}, \code{\link[rxode2:stat_cens]{geom_cens}}, \code{\link[rxode2]{ini}}, \code{\link[rxode2]{logit}}, \code{\link[rxode2:reexports]{lotri}}, \code{\link[rxode2]{model}}, \code{\link[rxode2]{probit}}, \code{\link[rxode2:probit]{probitInv}}, \code{\link[rxode2]{rxCat}}, \code{\link[rxode2]{rxClean}}, \code{\link[rxode2:rxSolve]{rxControl}}, \code{\link[rxode2:rxInits]{rxInit}}, \code{\link[rxode2]{rxLhs}}, \code{\link[rxode2]{rxModelVars}}, \code{\link[rxode2:rxModelVars]{rxModelVarsS3}}, \code{\link[rxode2:rxode2]{rxode}}, \code{\link[rxode2:rxode2]{RxODE}}, \code{\link[rxode2]{rxode2}}, \code{\link[rxode2:rxParams]{rxParam}}, \code{\link[rxode2]{rxParams}}, \code{\link[rxode2]{rxParams}}, \code{\link[rxode2]{rxSolve}}, \code{\link[rxode2]{rxSolve}}, \code{\link[rxode2]{rxState}}, \code{\link[rxode2]{stat_amt}}, \code{\link[rxode2]{stat_cens}}}
+  \item{rxode2}{\code{\link[rxode2:reexports]{add.dosing}}, \code{\link[rxode2:reexports]{add.sampling}}, \code{\link[rxode2]{as.rxUi}}, \code{\link[rxode2:reexports]{et}}, \code{\link[rxode2:reexports]{et}}, \code{\link[rxode2:reexports]{eventTable}}, \code{\link[rxode2:logit]{expit}}, \code{\link[rxode2:stat_amt]{geom_amt}}, \code{\link[rxode2:stat_cens]{geom_cens}}, \code{\link[rxode2]{ini}}, \code{\link[rxode2]{logit}}, \code{\link[rxode2:reexports]{lotri}}, \code{\link[rxode2]{model}}, \code{\link[rxode2]{probit}}, \code{\link[rxode2:probit]{probitInv}}, \code{\link[rxode2]{rxCat}}, \code{\link[rxode2]{rxClean}}, \code{\link[rxode2:rxSolve]{rxControl}}, \code{\link[rxode2:rxInits]{rxInit}}, \code{\link[rxode2]{rxLhs}}, \code{\link[rxode2]{rxModelVars}}, \code{\link[rxode2:rxModelVars]{rxModelVarsS3}}, \code{\link[rxode2:rxode2]{rxode}}, \code{\link[rxode2:rxode2]{RxODE}}, \code{\link[rxode2]{rxode2}}, \code{\link[rxode2:rxParams]{rxParam}}, \code{\link[rxode2]{rxParams}}, \code{\link[rxode2]{rxParams}}, \code{\link[rxode2]{rxSolve}}, \code{\link[rxode2]{rxSolve}}, \code{\link[rxode2]{rxState}}, \code{\link[rxode2]{stat_amt}}, \code{\link[rxode2]{stat_cens}}}
 }}
 
diff --git a/tests/testthat/test-nlm.R b/tests/testthat/test-nlm.R
new file mode 100644
index 00000000..f3fe06f7
--- /dev/null
+++ b/tests/testthat/test-nlm.R
@@ -0,0 +1,52 @@
+nmTest({
+  test_that("nlm makes sense", {
+
+
+    library(dplyr)
+    dsn <- data.frame(i=1:1000) %>%
+      mutate(time=exp(rnorm(n())),DV=rbinom(n(),1,exp(-1+time)/(1+exp(-1+time)))) %>%
+      mutate(id=1)
+
+    mod <- function() {
+      ini({
+        E0 <- 0.5
+        Em <- 0.5
+        E50 <- 2
+        g <- fix(2)
+      })
+      model({
+        v <- E0+Em*time^g/(E50^g+time^g)
+        ll(bin) ~ DV * v - log(1 + exp(v))
+      })
+    }
+
+    fit2 <- nlmixr(mod, dsn, est="nlm")
+
+    expect_true(inherits(fit2, "nlmixr2.nlm"))
+
+    fit3 <- fit2 %>% ini(g=unfix) %>% nlmixr2(dsn, "nlm", nlmControl(covMethod="nlm"))
+
+    expect_true(inherits(fit3, "nlmixr2.nlm"))
+
+    one.cmt <- function() {
+      ini({
+        tka <- 0.45
+        tcl <- log(c(0, 2.7, 100))
+        tv <- 3.45
+        add.sd <- 0.7
+      })
+      model({
+        ka <- exp(tka)
+        cl <- exp(tcl)
+        v <- exp(tv)
+        linCmt() ~ add(add.sd)
+      })
+    }
+
+    fit1 <- nlmixr(one.cmt, nlmixr2data::theo_sd, est="nlm")
+
+    expect_true(inherits(fit1, "nlmixr2.nlm"))
+
+  })
+
+})