diff --git a/NAMESPACE b/NAMESPACE
index 638bab91..23431ce2 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -336,10 +336,6 @@ export(compute_EIR_full)
 export(compute_NI)
 export(compute_NI_ix)
 export(compute_beta)
-export(compute_fqM)
-export(compute_fqM_ix)
-export(compute_fqZ)
-export(compute_fqZ_ix)
 export(compute_kappa)
 export(compute_local_frac)
 export(compute_terms)
diff --git a/R/blood_feeding.R b/R/blood_feeding.R
index 60f613f4..36c17372 100644
--- a/R/blood_feeding.R
+++ b/R/blood_feeding.R
@@ -67,9 +67,10 @@ setup_BFpar_static <- function(pars){
 #' @return none
 #' @export
 setup_BloodFeeding <- function(pars, i, s=1, BFopts = list(), residence=1, searchWts=1, F_circadian=NULL){
-  pars$BFpar$searchWts[[i]][[s]] = checkIt(searchWts, pars$Hpar[[i]]$nStrata)
-  pars$BFpar$relativeBitingRate[[i]][[s]] = checkIt(searchWts, pars$Hpar[[i]]$nStrata)
-  pars$BFpar$residence[[i]] = checkIt(residence, pars$Hpar[[i]]$nStrata)
+  nStrata = pars$Hpar[[i]]$nStrata
+  pars$BFpar$searchWts[[i]][[s]] = checkIt(searchWts, nStrata)
+  pars$BFpar$relativeBitingRate[[i]][[s]] = checkIt(searchWts, nStrata)
+  pars$BFpar$residence[[i]] = checkIt(residence, nStrata)
   if(is.null(F_circadian)) pars$BFpar$F_circadian[[s]] = function(d){return(1)}
   return(pars)
 }
diff --git a/R/compute_terms.R b/R/compute_terms.R
index a19ab1dd..50656dc7 100644
--- a/R/compute_terms.R
+++ b/R/compute_terms.R
@@ -82,15 +82,35 @@ compute_terms.cohort <- function(varslist, deout, pars, s, i) {
 #' @export
 compute_terms.human<- function(varslist, deout, pars, s, i) {
   time = deout[,1]
-  eir = c()
+
+  eir = list()
+  if(pars$nHosts>0) for(i in 1:pars$nHosts)
+    eir[[i]] = matrix(0, nrow = length(time), ncol = pars$Hpar[[i]]$nStrata)
+
+  fqZ = list()
+  for(s in 1:pars$nVectors)
+    fqZ[[s]] = matrix(0, nrow = length(time), ncol = pars$nPatches)
+
   for (ix in 1:length(time)){
-    yt = deout[i,-1]
-    pars = Transmission(time[i], yt, pars)
-    eir = c(eir, pars$EIR[[1]])
+    yt = deout[ix,-1]
+    pars = compute_vars_full(time[ix], yt, pars)
+
+    for(i in 1:pars$nHosts)
+      eir[[i]][ix,] = pars$EIR[[i]]
+
+    for(s in 1:pars$nVectors){
+      kappa[[s]][ix,] = pars$kappa[[i]]
+      fqZ[[s]][ix,] = F_fqZ(time[ix], yt, pars, s)
+    }
   }
-  ni = compute_NI(deout, pars)
-  fqZ = compute_fqZ(deout, pars)
-  pr = F_pr(varslist, pars)
+
+  pr = list()
+  ni = list()
+  for(i in 1:pars$nHosts){
+    ni[[i]] = compute_NI(deout, pars, i)
+    pr[[i]] = F_pr(varslist, pars, i)
+  }
+
   return(list(time=time,eir=eir,pr=pr,ni=ni,fqZ=fqZ))
 }
 
@@ -119,74 +139,11 @@ compute_terms_steady<- function(varslist, y_eq, pars, s, i) {
   pars = Transmission(0, y_eq, pars)
   eir = pars$EIR[[1]]
   fqZ <- F_fqZ(0, y_eq, pars, s)
-  ni <- F_X(0, y_eq, pars, i)/varslist$XH$H
+  ni <- F_X(0, y_eq, pars, i)/varslist$XH[[i]]$H
   pr <- F_pr(varslist, pars, i)
   return(list(eir=eir,pr=pr,kappa=kappa,fqZ=fqZ,ni=ni))
 }
 
-#' @title Compute the fqZ
-#' @description Using the output of [deSolve::ode] or [deSolve::dede],
-#' compute the fqZ at every point in time
-#' @param deout a matrix, the output of deSolve
-#' @param pars a [list]
-#' @param s the vector species index
-#' @return [vector]
-#' @export
-compute_fqZ <- function(deout, pars, s) {
-  d1 = length(deout[,1])
-  fqZ = sapply(1:d1, compute_fqZ_ix, deout=deout, pars=pars, s=s)
-  fqZ = shapeIt(fqZ, d1, pars$nPatches)
-  return(fqZ)
-}
-
-#' @title Compute the fqZ for the ith
-#' @description Using the output of [deSolve::ode] or [deSolve::dede],
-#' compute the fqZ at one point in time
-#' @param ix an [integer]
-#' @param deout a matrix, the output of deSolve
-#' @param pars a [list]
-#' @param s the vector species index
-#' @return [numeric]
-#' @export
-compute_fqZ_ix <- function(ix, deout, pars, s) {
-  t = deout[ix,1]
-  y = deout[ix,-1]
-  fqZ = F_fqZ(t, y, pars, s)
-  return(fqZ)
-}
-
-#' @title Compute the fqM
-#' @description Using the output of [deSolve::ode] or [deSolve::dede],
-#' compute the fqM at every point in time
-#' @param deout a matrix, the output of deSolve
-#' @param pars a [list]
-#' @param s the vector species index
-#' @return [vector]
-#' @export
-compute_fqM <- function(deout, pars, s) {
-  d1 = length(deout[,1])
-  fqM = sapply(1:d1, compute_fqM_ix, deout=deout, pars=pars, s=s)
-  fqM = shapeIt(fqM, d1, pars$nPatches)
-  return(fqM)
-}
-
-#' @title Compute the fqM for the ith
-#' @description Using the output of [deSolve::ode] or [deSolve::dede],
-#' compute the fqM at one point in time
-#' @param ix an [integer]
-#' @param deout a matrix, the output of deSolve
-#' @param pars a [list]
-#' @param s the vector species index
-#' @return [numeric]
-#' @export
-compute_fqM_ix <- function(ix, deout, pars, s) {
-  t = deout[ix,1]
-  y = deout[ix,-1]
-  fqM = F_fqM(t, y, pars, s)
-  return(fqM)
-}
-
-
 #' @title Compute the NI
 #' @description Using the output of [deSolve::ode] or [deSolve::dede],
 #' compute the NI for each stratum
@@ -198,7 +155,7 @@ compute_fqM_ix <- function(ix, deout, pars, s) {
 compute_NI <- function(deout, pars, i) {
   d1 = length(deout[,1])
   NI = sapply(1:d1, compute_NI_ix, deout=deout, pars=pars, i=i)
-  NI = shapeIt(NI, d1, pars$nStrata)
+  NI = shapeIt(NI, d1, pars$Hpar[[i]]$nStrata)
   return(NI)
 }
 
diff --git a/R/time_spent.R b/R/time_spent.R
index 184147f5..2a58e5dd 100644
--- a/R/time_spent.R
+++ b/R/time_spent.R
@@ -19,7 +19,7 @@ make_TimeSpent = function(pars, i, TimeSpent, opts = list()){
 #' @export
 make_TimeSpent.athome = function(pars, i, TimeSpent = "athome", opts = list()){
 
-    residence = pars$BFpar$residence[[i]]
+  residence = pars$BFpar$residence[[i]]
   TiSp = make_TimeSpent_athome(pars$nPatches, residence, opts)
   pars$BFpar$TimeSpent[[i]] = TiSp
   for(s in 1:pars$nVectors) pars = make_TaR(t, pars, i, s)
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 2ac0767c..2e6ea213 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -809,10 +809,6 @@ reference:
   - compute_terms.human
   - compute_terms.na
   - compute_terms_steady
-  - compute_fqZ
-  - compute_fqZ_ix
-  - compute_fqM
-  - compute_fqM_ix
   - compute_NI
   - compute_NI_ix
   - compute_vars_full
diff --git a/man/compute_fqM.Rd b/man/compute_fqM.Rd
deleted file mode 100644
index 3b060914..00000000
--- a/man/compute_fqM.Rd
+++ /dev/null
@@ -1,22 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/compute_terms.R
-\name{compute_fqM}
-\alias{compute_fqM}
-\title{Compute the fqM}
-\usage{
-compute_fqM(deout, pars, s)
-}
-\arguments{
-\item{deout}{a matrix, the output of deSolve}
-
-\item{pars}{a \link{list}}
-
-\item{s}{the vector species index}
-}
-\value{
-\link{vector}
-}
-\description{
-Using the output of \link[deSolve:ode]{deSolve::ode} or \link[deSolve:dede]{deSolve::dede},
-compute the fqM at every point in time
-}
diff --git a/man/compute_fqM_ix.Rd b/man/compute_fqM_ix.Rd
deleted file mode 100644
index 2c8a7018..00000000
--- a/man/compute_fqM_ix.Rd
+++ /dev/null
@@ -1,24 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/compute_terms.R
-\name{compute_fqM_ix}
-\alias{compute_fqM_ix}
-\title{Compute the fqM for the ith}
-\usage{
-compute_fqM_ix(ix, deout, pars, s)
-}
-\arguments{
-\item{ix}{an \link{integer}}
-
-\item{deout}{a matrix, the output of deSolve}
-
-\item{pars}{a \link{list}}
-
-\item{s}{the vector species index}
-}
-\value{
-\link{numeric}
-}
-\description{
-Using the output of \link[deSolve:ode]{deSolve::ode} or \link[deSolve:dede]{deSolve::dede},
-compute the fqM at one point in time
-}
diff --git a/man/compute_fqZ.Rd b/man/compute_fqZ.Rd
deleted file mode 100644
index f8e3a992..00000000
--- a/man/compute_fqZ.Rd
+++ /dev/null
@@ -1,22 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/compute_terms.R
-\name{compute_fqZ}
-\alias{compute_fqZ}
-\title{Compute the fqZ}
-\usage{
-compute_fqZ(deout, pars, s)
-}
-\arguments{
-\item{deout}{a matrix, the output of deSolve}
-
-\item{pars}{a \link{list}}
-
-\item{s}{the vector species index}
-}
-\value{
-\link{vector}
-}
-\description{
-Using the output of \link[deSolve:ode]{deSolve::ode} or \link[deSolve:dede]{deSolve::dede},
-compute the fqZ at every point in time
-}
diff --git a/man/compute_fqZ_ix.Rd b/man/compute_fqZ_ix.Rd
deleted file mode 100644
index 31ee274a..00000000
--- a/man/compute_fqZ_ix.Rd
+++ /dev/null
@@ -1,24 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/compute_terms.R
-\name{compute_fqZ_ix}
-\alias{compute_fqZ_ix}
-\title{Compute the fqZ for the ith}
-\usage{
-compute_fqZ_ix(ix, deout, pars, s)
-}
-\arguments{
-\item{ix}{an \link{integer}}
-
-\item{deout}{a matrix, the output of deSolve}
-
-\item{pars}{a \link{list}}
-
-\item{s}{the vector species index}
-}
-\value{
-\link{numeric}
-}
-\description{
-Using the output of \link[deSolve:ode]{deSolve::ode} or \link[deSolve:dede]{deSolve::dede},
-compute the fqZ at one point in time
-}