diff --git a/R/kron_tmp.R b/R/kron_tmp.R
deleted file mode 100644
index 49afce8f..00000000
--- a/R/kron_tmp.R
+++ /dev/null
@@ -1,8 +0,0 @@
-## Function translated using 'matlab.to.r()'
-## Then manually adjusted to make work
-## Author: Andrew Hooker
-
-kron_tmp <- function(mat1,mat2){
-    mat <- mat1 %x% mat2
-    return(mat)
-}
diff --git a/R/m4.R b/R/m4.R
deleted file mode 100644
index af727562..00000000
--- a/R/m4.R
+++ /dev/null
@@ -1,17 +0,0 @@
-## Function translated automatically using 'matlab.to.r()'
-## Author: Andrew Hooker
-
-m4 <- function(mv,n){
-#
-# size: (samps per individ^2 x samps per individ^2)
-  ns=n^2
-  mm4=zeros(ns,ns)
-  for(ct1 in 1:n){
-    for(ct2 in 1:n){
-      mm4[((ct1-1)*n+1):((ct1-1)*n+n),((ct2-1)*n+1):((ct2-1)*n+n)]= mv[ct1,ct2]*mv*2
-   }
-  }
-  ret=mm4
-return( ret ) 
-}
-
diff --git a/R/mf.R b/R/mf.R
deleted file mode 100644
index 149339ac..00000000
--- a/R/mf.R
+++ /dev/null
@@ -1,107 +0,0 @@
-#' The full  Fisher Information Matrix (FIM) for one individual
-#' 
-#' Compute the full  FIM for one individual given specific model(s), parameters, design and methods. 
-#' This computation makes no assumption that fixed and random effects are uncorrelated.  
-#' 
-#' @param xt_ind A vector of sample times.  
-#' @param model_switch A vector that is the same size as xt, specifying which model each sample belongs to.
-#' @param x A vector for the discrete design variables.
-#' @param a A vector of covariates.  
-#' 
-#' @inheritParams mftot
-#' 
-#' @return As a list:
-#' \item{ret}{The FIM for one individual}
-#' \item{poped.db}{A PopED database}
-#' 
-#' @seealso Used by \code{\link{mftot0}}.  
-#'  
-#' @family FIM
-#' @example tests/testthat/examples_fcn_doc/warfarin_basic.R
-#' @example tests/testthat/examples_fcn_doc/examples_mf.R
-#' 
-#' @export
-#' @keywords internal
-## Function translated automatically using 'matlab.to.r()'
-## Author: Andrew Hooker
-
-mf <- function(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db){
-  
-  numnotfixed_bpop = sum(poped.db$parameters$notfixed_bpop)
-  numnotfixed_d    = sum(poped.db$parameters$notfixed_d)
-  numnotfixed_covd = sum(poped.db$parameters$notfixed_covd)
-  numnotfixed_docc  = sum(poped.db$parameters$notfixed_docc)
-  numnotfixed_covdocc  = sum(poped.db$parameters$notfixed_covdocc)
-  numnotfixed_sigma  = sum(poped.db$parameters$notfixed_sigma)
-  numnotfixed_covsigma  = sum(poped.db$parameters$notfixed_covsigma)
-  
-  n=size(xt_ind,1)
-  ret = 0
-  
-  for(i in 1:poped.db$settings$iFOCENumInd){
-    b_ind = poped.db$parameters$b_global[,i,drop=F]
-    bocc_ind = poped.db$parameters$bocc_global[[i]]
-    
-    if((poped.db$settings$bCalculateEBE) ){#Calculate an EBE
-      epsi0 = zeros(1,length(poped.db$parameters$notfixed_sigma))
-      g=feval(poped.db$model$fg_pointer,x,a,bpop,b_ind,bocc_ind)
-      returnArgs <- feval(poped.db$model$ferror_pointer,model_switch,xt_ind,g,epsi0,poped.db) 
-      mean_data <- returnArgs[[1]]
-      poped.db <- returnArgs[[2]]
-      start_bind = t(b_ind)
-      b_ind = ind_estimates(mean_data,bpop,d,sigma,start_bind,(poped.db$settings$iApproximationMethod==2),FALSE,model_switch,xt_ind,x,a,b_ind,bocc_ind,poped.db)
-      #        b_ind2 = ind_estimates(mean_data,bpop,d,sigma,t(b_ind),(poped.db$settings$iApproximationMethod==2),FALSE,model_switch,xt_ind,x,a,b_ind,bocc_ind,poped.db)
-      
-      #b_ind2 = ind_estimates(mean_data,bpop,d,sigma,t(zeros(size(b_ind)[1],size(b_ind)[2])),!(poped.db$settings$iApproximationMethod==2),FALSE,model_switch,xt_ind,x,a,b_ind,bocc_ind,poped.db)
-      poped.db$mean_data = mean_data
-    }
-    
-    f1=zeros(n+n*n,numnotfixed_bpop+numnotfixed_d+numnotfixed_covd+numnotfixed_docc+numnotfixed_covdocc+numnotfixed_sigma+numnotfixed_covsigma)
-    returnArgs <- m1(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,d,poped.db) 
-    m1_tmp <- returnArgs[[1]]
-    poped.db <- returnArgs[[2]]
-    f1[1:n,1:numnotfixed_bpop]=m1_tmp
-    returnArgs <- m2(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,poped.db) 
-    f1[(n+1):(n+n*n),1:numnotfixed_bpop] <- returnArgs[[1]]
-    poped.db <- returnArgs[[2]]
-    if((numnotfixed_d+numnotfixed_covd+numnotfixed_docc+numnotfixed_covdocc+numnotfixed_sigma+numnotfixed_covsigma>0)){
-      returnArgs <- m3(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,TRUE,poped.db) 
-      f1[(n+1):(n+n*n),(numnotfixed_bpop+1):(numnotfixed_bpop+numnotfixed_d+numnotfixed_covd+numnotfixed_docc+numnotfixed_covdocc+numnotfixed_sigma+numnotfixed_covsigma)] <- returnArgs[[1]]
-      poped.db <- returnArgs[[2]]
-    }
-    f2=zeros(n+n*n,n+n*n)
-    returnArgs <-  v(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,poped.db) 
-    v_tmp <- returnArgs[[1]]
-    poped.db <- returnArgs[[2]]
-    if(any(v_tmp!=0)){
-      
-      v_tmp_inv = inv(v_tmp,pseudo_on_fail = T)
-      f2[1:n,1:n] = v_tmp_inv
-      
-      #tmp_m4_inv=0.25*m4(v_tmp_inv,n)
-      tmp_m4_inv <- 1/2*kronecker(v_tmp_inv,v_tmp_inv)
-      f2[(n+1):(n+n*n),(n+1):(n+n*n)] = tmp_m4_inv
-      
-      #f2[1:n,1:n]=v_tmp\diag_matlab(1,size(v_tmp))
-      #v_tmp\diag_matlab(1,size(v_tmp))
-      #f2[1:n,1:n]=inv(v_tmp)
-      #m4_tmp = m4(v_tmp,n)
-      #f2[(n+1):(n+n*n),(n+1):(n+n*n)]=m4_tmp\diag_matlab(1,size(m4_tmp))
-      #f2[(n+1):(n+n*n),(n+1):(n+n*n)]=inv(m4(v_tmp,n))
-    }
-    
-    if((sum(sum(f2!=0))==0) ){#Only zeros in f2, uses FIM = m1'*m1
-      ret = ret+t(m1_tmp)%*%m1_tmp
-    } else {
-      ret=ret+t(f1)%*%f2%*%f1
-    }
-    
-  }
-  #ret
-  ret = ret/poped.db$settings$iFOCENumInd
-  
-  
-  
-  return(list( ret= ret,poped.db=poped.db)) 
-}
-
diff --git a/R/mf3.R b/R/mf3.R
index c2d58f17..54d25510 100644
--- a/R/mf3.R
+++ b/R/mf3.R
@@ -32,6 +32,9 @@ mf3 <- function(model_switch,xt,x,a,bpop,d,sigma,docc,poped.db){
   numnotfixed_sigma  = sum(poped.db$parameters$notfixed_sigma)
   numnotfixed_covsigma  = sum(poped.db$parameters$notfixed_covsigma)
   
+  n_fixed_eff <- numnotfixed_bpop
+  n_rand_eff <- numnotfixed_d+numnotfixed_covd+numnotfixed_docc+numnotfixed_covdocc+numnotfixed_sigma+numnotfixed_covsigma
+
   n=size(xt,1)
   ret = 0
   
@@ -53,43 +56,79 @@ mf3 <- function(model_switch,xt,x,a,bpop,d,sigma,docc,poped.db){
       poped.db$mean_data = mean_data
     }
     
-    n_fixed_eff <- numnotfixed_bpop
-    n_rand_eff <- numnotfixed_d+numnotfixed_covd+numnotfixed_docc+numnotfixed_covdocc+numnotfixed_sigma+numnotfixed_covsigma
-    start_col <- 1
-    
-    f1=zeros(n+n*n,n_fixed_eff+n_rand_eff)
-    
     if(n_fixed_eff!=0){
-      returnArgs <- m1(model_switch,xt,x,a,bpop,b_ind,bocc_ind,d,poped.db) 
-      f1[1:n,start_col:(start_col+n_fixed_eff-1)] <- returnArgs[[1]]
+      returnArgs <- m1(model_switch,xt,x,a,bpop,b_ind,bocc_ind,d,poped.db)
+      m1_tmp <- returnArgs[[1]]
       poped.db <- returnArgs[[2]]
-      start_col <- start_col + n_fixed_eff
     }
     if(n_rand_eff!=0){
-      returnArgs <- m3(model_switch,xt,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,TRUE,poped.db) 
-      f1[(n+1):(n+n*n),start_col:(start_col+n_rand_eff-1)] <- returnArgs[[1]]
+      bUseVarSigmaDerivative = poped.db$settings$iFIMCalculationType != 4
+      returnArgs <- m3(model_switch,xt,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,bUseVarSigmaDerivative,poped.db)
+      m3_tmp <- returnArgs[[1]]
       poped.db <- returnArgs[[2]]
     }
-    
-    f2=zeros(n+n*n,n+n*n)
+    if(n_fixed_eff!=0 & n_rand_eff!=0 & poped.db$settings$iFIMCalculationType %in% c(0,5,6)){
+      returnArgs <- m2(model_switch,xt,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,poped.db) 
+      m2_tmp <- returnArgs[[1]]
+      poped.db <- returnArgs[[2]]
+    } else {
+      m2_tmp <- 0
+    }
+
     returnArgs <-  v(model_switch,xt,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,poped.db) 
     v_tmp <- returnArgs[[1]]
     poped.db <- returnArgs[[2]]
-    if(any(v_tmp!=0)){#If there are some non-zero elements to v_tmp
-      v_tmp_inv = inv(v_tmp,pseudo_on_fail = T)
-      f2[1:n,1:n] = v_tmp_inv
+
+    if (!poped.db$settings$iFIMCalculationType %in% c(5,7)) {
+      f1=zeros(n+n*n,n_fixed_eff+n_rand_eff)
+      f1[1:n,1:n_fixed_eff] <- m1_tmp
+      f1[(n+1):(n+n*n),1:n_fixed_eff] <- m2_tmp
+      f1[(n+1):(n+n*n),(n_fixed_eff+1):(n_fixed_eff+n_rand_eff)] <- m3_tmp
       
-      #tmp_m4_inv=0.25*m4(v_tmp_inv,n)
-      tmp_m4_inv <- 1/2*kronecker(v_tmp_inv,v_tmp_inv)
-      f2[(n+1):(n+n*n),(n+1):(n+n*n)] = tmp_m4_inv
-    }
-    if((all(f2==0))){
-      ret = ret+t(f1)%*%f1
+      if(any(v_tmp!=0)){#If there are some non-zero elements to v_tmp
+        f2=zeros(n+n*n,n+n*n)
+        
+        v_tmp_inv = inv(v_tmp,pseudo_on_fail = T)
+        f2[1:n,1:n] = v_tmp_inv
+        
+        tmp_m4_inv <- 1/2*kronecker(v_tmp_inv,v_tmp_inv)
+        f2[(n + 1):(n + n*n), (n + 1):(n + n*n)] = tmp_m4_inv
+        ret = ret + t(f1) %*% f2 %*% f1
+      } else {
+        ret = ret + t(f1) %*% f1
+      }
     } else {
-      ret = ret+t(f1)%*%f2%*%f1
+      v_tmp_inv = inv(v_tmp)
+      dim(m3_tmp) = c(n,n,n_rand_eff)
+      
+      tmp_fim = zeros(n_fixed_eff + n_rand_eff, n_fixed_eff + n_rand_eff)
+      tmp_fim[1:n_fixed_eff,1:n_fixed_eff] = 2*t(m1_tmp) %*% v_tmp_inv %*% m1_tmp
+
+      if (is.matrix(m2_tmp)) {
+        dim(m2_tmp) = c(n,n,n_fixed_eff)
+        for (m in 1:n_fixed_eff) {
+          for (k in 1:n_fixed_eff) {
+            tmp_fim[m,k] = tmp_fim[m,k] + trace_matrix(m2_tmp[,,m]%*%v_tmp_inv%*%m2_tmp[,,k]%*%v_tmp_inv)
+          }
+        }
+        for(m in 1:n_rand_eff){
+          for(k in 1:n_fixed_eff){
+            num = trace_matrix(m3_tmp[,,m]%*%v_tmp_inv%*%m2_tmp[,,k]%*%v_tmp_inv)
+            tmp_fim[n_fixed_eff + m, k]=num
+            tmp_fim[k, n_fixed_eff + m]=num
+          }
+        }
+      }
+      
+      for (m in 1:n_rand_eff) {
+        for (k in 1:n_rand_eff) {
+          tmp_fim[n_fixed_eff + m, n_fixed_eff + k] = trace_matrix(m3_tmp[,,m] %*% v_tmp_inv %*% m3_tmp[,,k] %*% v_tmp_inv)
+        }
+      }
+      ret = ret + 1/2*tmp_fim
     }
   }
   ret = ret/poped.db$settings$iFOCENumInd
   
-  return(list( ret= ret,poped.db=poped.db)) 
+  return(list(ret = ret,poped.db = poped.db)) 
 }
diff --git a/R/mf5.R b/R/mf5.R
deleted file mode 100644
index 31a21785..00000000
--- a/R/mf5.R
+++ /dev/null
@@ -1,83 +0,0 @@
-#' The reduced Fisher Information Matrix (FIM) for one individual, using the SD of RUV as a parameter.
-#'  
-#' 
-#' Compute the reduced FIM for one individual using the standard deviation of the residual unexplained variability (RUV) terms as a parameter, 
-#' given specific model(s), parameters, design and methods. 
-#' This computation 
-#' assumes that there is no correlation in the FIM between the fixed and random effects, 
-#' and set these elements in the FIM to zero.
-#' In addition all derivatives in the computation are made 
-#' with respect to the standard deviation of the RUV terms (sqrt(SIGMA) in NONMEM). 
-#' This matches what is done in PFIM, and assumes that the standard deviation of the residual unexplained variation is the estimated parameter
-#' (NOTE: NONMEM estimates the variance of the resudual unexplained variation by default).
-#' 
-#' @inheritParams mf3
-#' 
-#' @return As a list:
-#' \item{ret}{The FIM for one individual}
-#' \item{poped.db}{A PopED database}
-#' 
-#' @seealso Used by \code{\link{mftot4}}.  
-#' @family FIM
-#' 
-#' @example tests/testthat/examples_fcn_doc/warfarin_basic.R
-#' @example tests/testthat/examples_fcn_doc/examples_mf5.R
-#' @export
-#' @keywords internal
-## Function translated automatically using 'matlab.to.r()'
-## Author: Andrew Hooker
-
-mf5 <- function(model_switch,xt,x,a,bpop,d,sigma,docc,poped.db){
-  #==== Reduced FIM with derivative of variance w$r.t. the SD, NOT the
-  # variance
-  
-  
-  numnotfixed_bpop = sum(poped.db$parameters$notfixed_bpop)
-  numnotfixed_d    = sum(poped.db$parameters$notfixed_d)
-  numnotfixed_covd = sum(poped.db$parameters$notfixed_covd)
-  numnotfixed_docc  = sum(poped.db$parameters$notfixed_docc)
-  numnotfixed_covdocc  = sum(poped.db$parameters$notfixed_covdocc)
-  numnotfixed_sigma  = sum(poped.db$parameters$notfixed_sigma)
-  numnotfixed_covsigma  = sum(poped.db$parameters$notfixed_covsigma)
-  
-  
-  n=size(xt,1)
-  ret = 0
-  
-  for(i in 1:poped.db$settings$iFOCENumInd){
-    b_ind = poped.db$parameters$b_global[,i,drop=F]
-    bocc_ind = poped.db$parameters$bocc_global[[i]]
-    f1=zeros(n+n*n,numnotfixed_bpop+numnotfixed_d+numnotfixed_covd+numnotfixed_docc+numnotfixed_covdocc+numnotfixed_sigma+numnotfixed_covsigma)
-    returnArgs <- m1(model_switch,xt,x,a,bpop,b_ind,bocc_ind,d,poped.db) 
-    f1[1:n,1:numnotfixed_bpop] <- returnArgs[[1]]
-    poped.db <- returnArgs[[2]]
-    #w$r.t. the variance
-    #[f1[(n+1):(n+n*n),(numnotfixed_bpop+1):(numnotfixed_bpop+numnotfixed_d+numnotfixed_covd+numnotfixed_docc+numnotfixed_covdocc+numnotfixed_sigma+numnotfixed_covsigma)],poped.db]=m3(model_switch,xt,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,TRUE,poped.db)
-    #w$r.t. the sd
-    returnArgs <- m3(model_switch,xt,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,FALSE,poped.db) 
-    f1[(n+1):(n+n*n),(numnotfixed_bpop+1):(numnotfixed_bpop+numnotfixed_d+numnotfixed_covd+numnotfixed_docc+numnotfixed_covdocc+numnotfixed_sigma+numnotfixed_covsigma)] <- returnArgs[[1]]
-    poped.db <- returnArgs[[2]]
-    
-    f2=zeros(n+n*n,n+n*n)
-    returnArgs <-  v(model_switch,xt,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,poped.db) 
-    v_tmp <- returnArgs[[1]]
-    poped.db <- returnArgs[[2]]
-    if(any(v_tmp!=0)){#If the inverse is not empty
-      
-      v_tmp_inv = inv(v_tmp,pseudo_on_fail = T)
-      f2[1:n,1:n] = v_tmp_inv
-      
-      #tmp_m4_inv=0.25*m4(v_tmp_inv,n)
-      tmp_m4_inv <- 1/2*kronecker(v_tmp_inv,v_tmp_inv)
-      f2[(n+1):(n+n*n),(n+1):(n+n*n)] = tmp_m4_inv
-      
-      #f2[1:n,1:n]=inv(v_tmp)
-      #tmp_m4=m4(v_tmp,n)
-      #f2[(n+1):(n+n*n),(n+1):(n+n*n)]=inv(tmp_m4)
-    }
-    ret = ret+t(f1)%*%f2%*%f1
-  }
-  ret = ret/poped.db$settings$iFOCENumInd
-  return(list( ret= ret,poped.db=poped.db)) 
-}
-
diff --git a/R/mf6.R b/R/mf6.R
deleted file mode 100644
index 54046b28..00000000
--- a/R/mf6.R
+++ /dev/null
@@ -1,107 +0,0 @@
-#' The full Fisher Information Matrix (FIM) for one individual parameterized with A,B,C matrices & using the derivative of variance.
-#' 
-#' Compute the full FIM for one individual given specific model(s), parameters, design and methods. 
-#' This computation parameterizes the FIM calculation using 
-#' A,B,C matrices (as in Retout \emph{et al.}) but uses the derivative of variances.
-#' Should give the same answer as \code{\link{mf}} but computation times may be different.   
-#' 
-#' @inheritParams mf
-#' 
-#' @return As a list:
-#' \item{ret}{The FIM for one individual}
-#' \item{poped.db}{A PopED database}
-#' 
-#' @seealso Used by \code{\link{mftot5}}.  
-#'  
-#' @family FIM
-#' @references S. Retout and F. Mentre, "Further developments of the Fisher Information Matrix in
-#' nonlinear mixed effects models with evaluation in population pharmacokinetics", J. of Biopharm. Stats., 13(2), 2003.
-#' 
-#' @example tests/testthat/examples_fcn_doc/warfarin_basic.R
-#' @example tests/testthat/examples_fcn_doc/examples_mf6.R
-#' @export
-#' @keywords internal
-## Function translated automatically using 'matlab.to.r()'
-## Author: Andrew Hooker
-
-mf6 <- function(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db){
-  
-  #Calculate FIM with another parameterization, i$e. the parametrization used
-  #in Retout et al but with derivative of variance instead of SD for sigma
-  
-  numnotfixed_bpop = sum(poped.db$parameters$notfixed_bpop)
-  numnotfixed_d    = sum(poped.db$parameters$notfixed_d)
-  numnotfixed_covd = sum(poped.db$parameters$notfixed_covd)
-  numnotfixed_docc  = sum(poped.db$parameters$notfixed_docc)
-  numnotfixed_covdocc  = sum(poped.db$parameters$notfixed_covdocc)
-  numnotfixed_sigma  = sum(poped.db$parameters$notfixed_sigma)
-  numnotfixed_covsigma  = sum(poped.db$parameters$notfixed_covsigma)
-  
-  poped.db$settings$bCalculateEBE = FALSE
-  
-  n=size(xt_ind,1)
-  ret = 0
-  
-  for(i in 1:poped.db$settings$iFOCENumInd){
-    b_ind = poped.db$parameters$b_global[,i,drop=F]
-    bocc_ind = poped.db$parameters$bocc_global[[i]]
-    
-    if((poped.db$settings$bCalculateEBE) ){#Calculate an EBE
-      epsi0 = zeros(1,length(poped.db$parameters$notfixed_sigma))
-      g=feval(poped.db$model$fg_pointer,x,a,bpop,b_ind,bocc_ind)
-      returnArgs <- feval(poped.db$model$ferror_pointer,model_switch,xt_ind,g,epsi0,poped.db) 
-      mean_data <- returnArgs[[1]]
-      poped.db <- returnArgs[[2]]
-      start_bind = t(b_ind)
-      b_ind = ind_estimates(mean_data,bpop,d,sigma,start_bind,(poped.db$settings$iApproximationMethod==2),FALSE,model_switch,xt_ind,x,a,b_ind,bocc_ind,poped.db)
-      #        b_ind2 = ind_estimates(mean_data,bpop,d,sigma,t(b_ind),(poped.db$settings$iApproximationMethod==2),FALSE,model_switch,xt_ind,x,a,b_ind,bocc_ind,poped.db)
-      
-      #b_ind2 = ind_estimates(mean_data,bpop,d,sigma,t(zeros(size(b_ind)[1],size(b_ind)[2])),!(poped.db$settings$iApproximationMethod==2),FALSE,model_switch,xt_ind,x,a,b_ind,bocc_ind,poped.db)
-      poped.db$mean_data = mean_data
-    }
-    
-    tmp_fim=zeros(numnotfixed_bpop+numnotfixed_d+numnotfixed_covd+numnotfixed_docc+numnotfixed_covdocc+numnotfixed_sigma+numnotfixed_covsigma,numnotfixed_bpop+numnotfixed_d+numnotfixed_covd+numnotfixed_docc+numnotfixed_covdocc+numnotfixed_sigma+numnotfixed_covsigma)
-    returnArgs <-  m1(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,d,poped.db) 
-    m1_tmp <- returnArgs[[1]]
-    poped.db <- returnArgs[[2]]
-    returnArgs <-  m2(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,poped.db) 
-    m2_tmp <- returnArgs[[1]]
-    poped.db <- returnArgs[[2]]
-    returnArgs <-  m3(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,TRUE,poped.db) 
-    m3_tmp <- returnArgs[[1]]
-    poped.db <- returnArgs[[2]]
-    returnArgs <-  v(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,poped.db) 
-    v_tmp <- returnArgs[[1]]
-    poped.db <- returnArgs[[2]]
-    invv = inv(v_tmp)
-    
-    numtotalnotfixed = numnotfixed_d+numnotfixed_covd+numnotfixed_docc+numnotfixed_covdocc+numnotfixed_sigma+numnotfixed_covsigma
-    dim(m2_tmp) = c(n,n,numnotfixed_bpop)
-    dim(m3_tmp) = c(n,n,numtotalnotfixed)
-
-    tmp_fim[1:numnotfixed_bpop,1:numnotfixed_bpop]=2*t(m1_tmp)%*%invv%*%m1_tmp
-    for(m in 1:numnotfixed_bpop){
-      for(k in 1:numnotfixed_bpop){
-        tmp_fim[m,k]=tmp_fim[m,k]+trace_matrix(m2_tmp[,,m]%*%invv%*%m2_tmp[,,k]%*%invv)
-      }
-    }
-    
-    for(m in 1:numtotalnotfixed){
-      for(k in 1:numtotalnotfixed){
-        tmp_fim[numnotfixed_bpop+m,numnotfixed_bpop+k]=trace_matrix(m3_tmp[,,m]%*%invv%*%m3_tmp[,,k]%*%invv)
-      }
-    }
-    
-    for(m in 1:numtotalnotfixed){
-      for(k in 1:numnotfixed_bpop){
-        num =trace_matrix(m3_tmp[,,m]%*%invv%*%m2_tmp[,,k]%*%invv)
-        tmp_fim[numnotfixed_bpop+m,k]=num
-        tmp_fim[k,numnotfixed_bpop+m]=num
-      }
-    }
-    ret = ret+1/2*tmp_fim
-  }
-  ret = ret/poped.db$settings$iFOCENumInd
-  return(list( ret= ret,poped.db=poped.db)) 
-}
-
diff --git a/R/mf8.R b/R/mf8.R
deleted file mode 100644
index 894d5959..00000000
--- a/R/mf8.R
+++ /dev/null
@@ -1,80 +0,0 @@
-#' The reduced Fisher Information Matrix (FIM) for one individual parameterized with A,B,C matrices & using the derivative of variance.
-#' 
-#' Compute the reduced FIM for one individual given specific model(s), parameters, design and methods. 
-#' This computation assumes that there is no correlation in the FIM between the fixed and random effects, 
-#' and set these elements in the FIM to zero.
-#' This computation parameterizes the FIM calculation using 
-#' A,B,C matrices (as in Retout \emph{et al.}) but uses the derivative of variances.
-#' Should give the same answer as \code{\link{mf3}} but computation times may be different.   
-#' 
-#' @inheritParams mf
-#' 
-#' @return As a list:
-#' \item{ret}{The FIM for one individual}
-#' \item{poped.db}{A PopED database}
-#' 
-#' @seealso Used by \code{\link{mftot7}}.  
-#' @family FIM
-#' 
-#' @references S. Retout and F. Mentre, "Further developments of the Fisher Information Matrix in
-#' nonlinear mixed effects models with evaluation in population pharmacokinetics", J. of Biopharm. Stats., 13(2), 2003.
-#' @example tests/testthat/examples_fcn_doc/warfarin_basic.R
-#' @example tests/testthat/examples_fcn_doc/examples_mf8.R
-#' @export
-#' @keywords internal
-## Function translated automatically using 'matlab.to.r()'
-## Author: Andrew Hooker
-
-mf8 <- function(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db){
-
-#Calculate FIM with another parameterization, i$e. the parametrization used
-#in Retout et al but with derivative of variance instead of SD for sigma and the reduced FIM
-
-numnotfixed_bpop = sum(poped.db$parameters$notfixed_bpop)
-numnotfixed_d    = sum(poped.db$parameters$notfixed_d)
-numnotfixed_covd = sum(poped.db$parameters$notfixed_covd)
-numnotfixed_docc  = sum(poped.db$parameters$notfixed_docc)
-numnotfixed_covdocc  = sum(poped.db$parameters$notfixed_covdocc)
-numnotfixed_sigma  = sum(poped.db$parameters$notfixed_sigma)
-numnotfixed_covsigma  = sum(poped.db$parameters$notfixed_covsigma)
-
-n=size(xt_ind,1)
-ret = 0
-
-for(i in 1:poped.db$settings$iFOCENumInd){
-    b_ind = poped.db$parameters$b_global[,i,drop=F]
-    bocc_ind = poped.db$parameters$bocc_global[[i]]
-    tmp_fim=zeros(numnotfixed_bpop+numnotfixed_d+numnotfixed_covd+numnotfixed_docc+numnotfixed_covdocc+numnotfixed_sigma+numnotfixed_covsigma,numnotfixed_bpop+numnotfixed_d+numnotfixed_covd+numnotfixed_docc+numnotfixed_covdocc+numnotfixed_sigma+numnotfixed_covsigma)
-     returnArgs <-  m1(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,d,poped.db) 
-m1_tmp <- returnArgs[[1]]
-poped.db <- returnArgs[[2]]
-     returnArgs <-  m3(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,TRUE,poped.db) 
-m3_tmp <- returnArgs[[1]]
-poped.db <- returnArgs[[2]]
-     returnArgs <-  v(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,poped.db) 
-v_tmp <- returnArgs[[1]]
-poped.db <- returnArgs[[2]]
-    invv = inv(v_tmp)
-    
-    numtotalnotfixed = numnotfixed_d+numnotfixed_covd+numnotfixed_docc+numnotfixed_covdocc+numnotfixed_sigma+numnotfixed_covsigma
-    dim(m3_tmp) = c(n,n,numtotalnotfixed)
-
-    tmp_fim[1:numnotfixed_bpop,1:numnotfixed_bpop]=2*t(m1_tmp)%*%invv%*%m1_tmp
-    for(m in 1:numnotfixed_bpop){
-        for(k in 1:numnotfixed_bpop){
-            tmp_fim[m,k]=tmp_fim[m,k]
-        }
-    }
-
-    for(m in 1:numtotalnotfixed){
-        for(k in 1:numtotalnotfixed){
-            tmp_fim[numnotfixed_bpop+m,numnotfixed_bpop+k]=trace_matrix(m3_tmp[,,m]%*%invv%*%m3_tmp[,,k]%*%invv)
-        }
-    }
-
-    ret = ret+1/2*tmp_fim
-}
-ret = ret/poped.db$settings$iFOCENumInd
-return(list( ret= ret,poped.db=poped.db)) 
-}
-
diff --git a/R/mf_all.R b/R/mf_all.R
index 57cedd53..369b6c7d 100644
--- a/R/mf_all.R
+++ b/R/mf_all.R
@@ -5,14 +5,14 @@
 mf_all <- function(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db){
 
     returnArgs <- switch(poped.db$settings$iFIMCalculationType+1,
-                         mf(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db), #Default (with no assumption that bpop and b are uncorrelated)
+                         mf3(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db), #Default (with no assumption that bpop and b are uncorrelated)
                          mf3(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db), #Reduced FIM
                          stop("Not yet implemented"), #Weighted models
                          stop("Not yet implemented"), #Loc models
-                         mf5(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db), #Reduced FIM with derivative of SD sigma
-                         mf6(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db), #FULL FIM parameterized with A,B,C matrices & derivative of variance
+                         mf3(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db), #Reduced FIM with derivative of SD sigma
+                         mf3(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db), #FULL FIM parameterized with A,B,C matrices & derivative of variance
                          mf7(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db), #Calculate one model switch at a time, good for large matrices
-                         mf8(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db)) #Reduced FIM parameterized with A,B,C matrices & derivative of variance
+                         mf3(model_switch,xt_ind,x,a,bpop,d,sigma,docc,poped.db)) #Reduced FIM parameterized with A,B,C matrices & derivative of variance
     
     if(is.null(returnArgs)) stop(sprintf('Unknown FIM-calculation type'))
     ret <- returnArgs[[1]]
diff --git a/R/mftot.R b/R/mftot.R
index fb8a05f5..8c34b8bc 100644
--- a/R/mftot.R
+++ b/R/mftot.R
@@ -30,19 +30,28 @@
 ## Author: Andrew Hooker
 
 mftot <- function(model_switch,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db){
-    
-    returnArgs <- switch(poped.db$settings$iFIMCalculationType+1,
-                         mftot0(model_switch,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db), 
-                         mftot1(model_switch,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db), 
-                         mftot2(model_switch,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db),
-                         mftot3(model_switch,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db), 
-                         mftot4(model_switch,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db),
-                         mftot5(model_switch,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db),
-                         mftot6(model_switch,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db),
-                         mftot7(model_switch,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db)) 
-    
-    if(is.null(returnArgs)) stop(sprintf('Unknown FIM-calculation type'))
-    ret <- returnArgs[[1]]
-    poped.db <- returnArgs[[2]]
-    return(list( ret= ret,poped.db =poped.db )) 
+  m=size(ni,1)
+  s=0
+  for(i in 1:m){
+    if((ni[i]!=0 && groupsize[i]!=0)){
+      if((!isempty(x))){
+        x_i = t(x[i,,drop=F])      
+      } else {
+        x_i =  zeros(0,1)
+      }
+      if((!isempty(a))){
+        a_i = t(a[i,,drop=F])
+      } else {
+        a_i =  zeros(0,1)
+      }
+      # mf_all <- function(model_switch,xt,x,a,bpop,d,sigma,docc,poped.db){
+      returnArgs <- mf_all(t(model_switch[i,1:ni[i,drop=F],drop=F]),t(xt[i,1:ni[i,drop=F],drop=F]),x_i,a_i,bpop,d,sigma,docc,poped.db) 
+      if(is.null(returnArgs)) stop(sprintf('Unknown FIM-calculation type'))
+      mf_tmp <- returnArgs[[1]]
+      poped.db <- returnArgs[[2]]
+      s=s+groupsize[i]*mf_tmp
+    }
+  }
+
+  return(list( ret= s,poped.db =poped.db )) 
 }
diff --git a/R/mftot0.R b/R/mftot0.R
deleted file mode 100644
index 6ddd71c5..00000000
--- a/R/mftot0.R
+++ /dev/null
@@ -1,46 +0,0 @@
-#' The full Fisher Information Matrix (FIM)
-#' 
-#' Compute the full FIM given specific model(s), parameters, design and methods. 
-#' This computation makes no assumption that fixed and random effects are uncorrelated.  
-#' 
-#' @inheritParams mftot
-#' 
-#' @return As a list:
-#' \item{ret}{The FIM}
-#' \item{poped.db}{A PopED database}
-#' 
-#' @seealso For an easier function to use, please see \code{\link{evaluate.fim}}.  
-#'  
-#' @family FIM
-#' @example tests/testthat/examples_fcn_doc/warfarin_optimize.R
-#' @example tests/testthat/examples_fcn_doc/examples_mftot0.R
-#' @export
-#' @keywords internal
-## Function translated automatically using 'matlab.to.r()'
-## Author: Andrew Hooker
-
-mftot0 <- function(model_switch,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db){
-m=size(ni,1)
-s=0
-for(i in 1:m){
-    if((ni[i]!=0 && groupsize[i]!=0)){
-       if((!isempty(x))){
-         x_i = t(x[i,,drop=F])      
-       } else {
-         x_i =  zeros(0,1)
-       }
-       if((!isempty(a))){
-        a_i = t(a[i,,drop=F])
-       } else {
-        a_i =  zeros(0,1)
-       }
-         returnArgs <- mf(t(model_switch[i,1:ni[i,drop=F],drop=F]),t(xt[i,1:ni[i,drop=F],drop=F]), 
-                            x_i,a_i,bpop,d,sigma,docc,poped.db) 
-mf_tmp <- returnArgs[[1]]
-poped.db <- returnArgs[[2]]
-        s=s+groupsize[i]*mf_tmp
-    }
-}
-ret = s
-return(list( ret= ret,poped.db =poped.db )) 
-}
diff --git a/R/mftot1.R b/R/mftot1.R
deleted file mode 100644
index 403e2fec..00000000
--- a/R/mftot1.R
+++ /dev/null
@@ -1,48 +0,0 @@
-#' The reduced Fisher Information Matrix (FIM)
-#' 
-#' Compute the reduced FIM given specific model(s), parameters, design and methods. 
-#' This computation assumes that there is no correlation in the FIM between the fixed and random effects, 
-#' and set these elements in the FIM to zero.
-#' 
-#' @inheritParams mftot
-#' 
-#' @return As a list:
-#' \item{ret}{The FIM}
-#' \item{poped.db}{A PopED database}
-#' 
-#' @seealso For an easier function to use, please see \code{\link{evaluate.fim}}.  
-#'  
-#' @family FIM
-#' 
-#' @example tests/testthat/examples_fcn_doc/warfarin_optimize.R
-#' @example tests/testthat/examples_fcn_doc/examples_mftot1.R
-#' @export
-#' @keywords internal
-## Function translated automatically using 'matlab.to.r()'
-## Author: Andrew Hooker
-
-mftot1 <- function(model_switch,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db){
-  m=size(ni,1)
-  s=0
-  for(i in 1:m){
-    if((ni[i]!=0 && groupsize[i]!=0)){
-      if((!isempty(x))){
-        x_i = t(x[i,,drop=F])      
-      } else {
-        x_i =  zeros(0,1)
-      }
-      if((!isempty(a))){
-        a_i = t(a[i,,drop=F])
-      } else {
-        a_i =  zeros(0,1)
-      }
-      returnArgs <- mf3(t(model_switch[i,1:ni[i,drop=F],drop=F]),t(xt[i,1:ni[i,drop=F],drop=F]),x_i,a_i,bpop,d,sigma,docc,poped.db) 
-      mf3_tmp <- returnArgs[[1]]
-      poped.db <- returnArgs[[2]]
-      s=s+groupsize[i]*mf3_tmp
-    }
-  }
-  
-  ret=s
-  return(list( ret= ret,poped.db =poped.db )) 
-}
diff --git a/R/mftot2.R b/R/mftot2.R
deleted file mode 100644
index 2a3a3534..00000000
--- a/R/mftot2.R
+++ /dev/null
@@ -1,21 +0,0 @@
-#' The Fisher Information Matrix (FIM) using weighted models
-#' 
-#' Compute the FIM using weighted models given specific model(s), parameters, design and methods. Not currently available.
-#' 
-#' @inheritParams mftot
-#' 
-#' @seealso For an easier function to use, please see \code{\link{evaluate.fim}}.  
-#'  
-#' @family FIM
-#' 
-#' @export
-#' @keywords internal
-
-## Function translated automatically using 'matlab.to.r()'
-## Author: Andrew Hooker
-
-mftot2 <- function(model_switch,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db){
-
-stop(sprintf('This iFIMCalculationType option in not currently used in PopED'))
-
-}
diff --git a/R/mftot3.R b/R/mftot3.R
deleted file mode 100644
index ca00ea36..00000000
--- a/R/mftot3.R
+++ /dev/null
@@ -1,22 +0,0 @@
-#' The Fisher Information Matrix (FIM) some other method
-#' 
-#' Compute the FIM using some other method given specific model(s), parameters, design and methods. This is a placeholder.
-#' 
-#' @inheritParams mftot
-#' 
-#' 
-#' @seealso For an easier function to use, please see \code{\link{evaluate.fim}}.  
-#'  
-#' @family FIM
-#' 
-#' @export
-#' @keywords internal
-
-## Function translated automatically using 'matlab.to.r()'
-## Author: Andrew Hooker
-
-mftot3 <- function(model_switch,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db){
-
-stop(sprintf('This iFIMCalculationType option in not currently used in PopED'))
-
-}
diff --git a/R/mftot4.R b/R/mftot4.R
deleted file mode 100644
index 9af842ee..00000000
--- a/R/mftot4.R
+++ /dev/null
@@ -1,58 +0,0 @@
-#' The reduced Fisher Information Matrix (FIM) using the SD of RUV as a parameter.
-#'  
-#' 
-#' Compute the reduced FIM using the standard deviation of the residual unexplained variability (RUV) terms as a parameter, 
-#' given specific model(s), parameters, design and methods. 
-#' This computation 
-#' assumes that there is no correlation in the FIM between the fixed and random effects, 
-#' and set these elements in the FIM to zero.
-#' In addition all derivatives in the computation are made 
-#' with respect to the standard deviation of the RUV terms (sqrt(SIGMA) in NONMEM). 
-#' This matches what is done in PFIM, and assumes that the standard deviation of the residual unexplained variation is the estimated parameter
-#' (NOTE: NONMEM estimates the variance of the resudual unexplained variation by default).
-#' 
-#' @inheritParams mftot
-#' 
-#' @return As a list:
-#' \item{ret}{The FIM}
-#' \item{poped.db}{A PopED database}
-#' 
-#' @seealso For an easier function to use, please see \code{\link{evaluate.fim}}.  
-#'  
-#' @family FIM
-#' 
-#' @example tests/testthat/examples_fcn_doc/warfarin_optimize.R
-#' @example tests/testthat/examples_fcn_doc/examples_mftot4.R
-#' @export
-#' @keywords internal
-## Function translated automatically using 'matlab.to.r()'
-## Author: Andrew Hooker
-
-## Function translated automatically using 'matlab.to.r()'
-## Author: Andrew Hooker
-
-mftot4 <- function(model_switch,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db){
-  m=size(ni,1)
-  s=0
-  for(i in 1:m){
-     if((ni[i]!=0 && groupsize[i]!=0)){
-          if((!isempty(x))){
-            x_i = t(x[i,,drop=F])      
-         } else {
-            x_i =  zeros(0,1)
-         }
-         if((!isempty(a))){
-            a_i = t(a[i,,drop=F])
-         } else {
-            a_i =  zeros(0,1)
-         }
-        
-       returnArgs <- mf5(t(model_switch[i,1:ni[i,drop=F],drop=F]),t(xt[i,1:ni[i,drop=F],drop=F]),x_i,a_i, bpop,d,sigma,docc,poped.db) 
-mf_tmp <- returnArgs[[1]]
-poped.db <- returnArgs[[2]]
-      s=s+groupsize[i]*mf_tmp
-    }
-  }
-  ret=s
-return(list( ret= ret,poped.db =poped.db )) 
-}
diff --git a/R/mftot5.R b/R/mftot5.R
deleted file mode 100644
index ad002931..00000000
--- a/R/mftot5.R
+++ /dev/null
@@ -1,52 +0,0 @@
-#' The full Fisher Information Matrix (FIM) parameterized with A,B,C matrices & using the derivative of variance.
-#' 
-#' Compute the full FIM given specific model(s), parameters, design and methods. 
-#' This computation parameterizes the FIM calculation using 
-#' A,B,C matrices (as in Retout \emph{et al.}) but uses the derivative of variances.
-#' Should give the same answer as \code{\link{mftot0}} but computation times may be different.   
-#' 
-#' @inheritParams mftot
-#' 
-#' @return As a list:
-#' \item{ret}{The FIM}
-#' \item{poped.db}{A PopED database}
-#' 
-#' @seealso For an easier function to use, please see \code{\link{evaluate.fim}}.  
-#'  
-#' @family FIM
-#' @references S. Retout and F. Mentre, "Further developments of the Fisher Information Matrix in
-#' nonlinear mixed effects models with evaluation in population pharmacokinetics", J. of Biopharm. Stats., 13(2), 2003.
-#' 
-#' @example tests/testthat/examples_fcn_doc/warfarin_optimize.R
-#' @example tests/testthat/examples_fcn_doc/examples_mftot5.R
-#' @export
-#' @keywords internal
-## Function translated automatically using 'matlab.to.r()'
-## Author: Andrew Hooker
-
-mftot5 <- function(model_switch,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db){
-m=size(ni,1)
-s=0
-for(i in 1:m){
-   if((ni[i]!=0 && groupsize[i]!=0)){
-      
-         if((!isempty(x))){
-            x_i = t(x[i,,drop=F])      
-         } else {
-            x_i =  zeros(0,1)
-         }
-         if((!isempty(a))){
-            a_i = t(a[i,,drop=F])
-         } else {
-            a_i =  zeros(0,1)
-         }
-
-         returnArgs <- mf6(t(model_switch[i,1:ni[i,drop=F],drop=F]),t(xt[i,1:ni[i,drop=F],drop=F]),x_i,a_i,bpop,d,sigma,docc,poped.db) 
-mf_tmp <- returnArgs[[1]]
-poped.db <- returnArgs[[2]]
-        s=s+groupsize[i]*mf_tmp
-    }
-}
-ret = s
-return(list( ret= ret,poped.db =poped.db )) 
-}
diff --git a/R/mftot6.R b/R/mftot6.R
deleted file mode 100644
index e01ebe31..00000000
--- a/R/mftot6.R
+++ /dev/null
@@ -1,48 +0,0 @@
-#' The full Fisher Information Matrix (FIM) Calculating one model switch at a time, good for large matrices.
-#' 
-#' Compute the full FIM given specific model(s), parameters, design and methods. 
-#' This computation calculates the FIM for each model switch separately.  Correlations between the models parameters are assumed to be zero.
-#' 
-#' @inheritParams mftot
-#' 
-#' @return As a list:
-#' \item{ret}{The FIM}
-#' \item{poped.db}{A PopED database}
-#' 
-#' @seealso For an easier function to use, please see \code{\link{evaluate.fim}}.  
-#' @family FIM
-#' 
-#' 
-#' @example tests/testthat/examples_fcn_doc/warfarin_optimize.R
-#' @example tests/testthat/examples_fcn_doc/examples_mftot6.R
-#' @export
-#' @keywords internal
-## Function translated automatically using 'matlab.to.r()'
-## Author: Andrew Hooker
-
-mftot6 <- function(model_switch,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db){
-m=size(ni,1)
-s=0
-for(i in 1:m){
-   if((ni[i]!=0 && groupsize[i]!=0)){
-      
-         if((!isempty(x))){
-           x_i = t(x[i,,drop=F])      
-         } else {
-           x_i =  zeros(0,1)
-         }
-         if((!isempty(a))){
-           a_i = t(a[i,,drop=F])
-         } else {
-           a_i =  zeros(0,1)
-         }
-      
-         returnArgs <- mf7(t(model_switch[i,1:ni[i,drop=F],drop=F]),t(xt[i,1:ni[i,drop=F],drop=F]),x_i,a_i,bpop,d,sigma,docc,poped.db) 
-mf_tmp <- returnArgs[[1]]
-poped.db <- returnArgs[[2]]
-        s=s+groupsize[i]*mf_tmp
-    }
-}
-ret = s
-return(list( ret= ret,poped.db =poped.db )) 
-}
diff --git a/R/mftot7.R b/R/mftot7.R
deleted file mode 100644
index c66c34ca..00000000
--- a/R/mftot7.R
+++ /dev/null
@@ -1,54 +0,0 @@
-#' The reduced Fisher Information Matrix (FIM) parameterized with A,B,C matrices & using the derivative of variance.
-#' 
-#' Compute the reduced FIM given specific model(s), parameters, design and methods. 
-#' This computation assumes that there is no correlation in the FIM between the fixed and random effects, 
-#' and set these elements in the FIM to zero.
-#' This computation parameterizes the FIM calculation using 
-#' A,B,C matrices (as in Retout \emph{et al.}) but uses the derivative of variances.
-#' Should give the same answer as \code{\link{mftot1}} but computation times may be different.   
-#' 
-#' @inheritParams mftot
-#' 
-#' @return As a list:
-#' \item{ret}{The FIM}
-#' \item{poped.db}{A PopED database}
-#' 
-#' @seealso For an easier function to use, please see \code{\link{evaluate.fim}}.  
-#' @family FIM
-#' 
-#' @references S. Retout and F. Mentre, "Further developments of the Fisher Information Matrix in
-#' nonlinear mixed effects models with evaluation in population pharmacokinetics", J. of Biopharm. Stats., 13(2), 2003.
-#' 
-#' @example tests/testthat/examples_fcn_doc/warfarin_optimize.R
-#' @example tests/testthat/examples_fcn_doc/examples_mftot7.R
-#' @export
-#' @keywords internal
-## Function translated automatically using 'matlab.to.r()'
-## Author: Andrew Hooker
-
-mftot7 <- function(model_switch,groupsize,ni,xt,x,a,bpop,d,sigma,docc,poped.db){
-m=size(ni,1)
-s=0
-for(i in 1:m){
-   if((ni[i]!=0 && groupsize[i]!=0)){
-
-         if((!isempty(x))){
-           x_i = t(x[i,,drop=F])
-         } else {
-           x_i =  zeros(0,1)
-         }
-         if((!isempty(a))){
-           a_i = t(a[i,,drop=F])
-         } else {
-           a_i =  zeros(0,1)
-         }
-
-         returnArgs <- mf8(t(model_switch[i,1:ni[i,drop=F],drop=F]),t(xt[i,1:ni[i,drop=F],drop=F]),x_i,a_i,bpop,d,sigma,docc,poped.db) 
-mf_tmp <- returnArgs[[1]]
-poped.db <- returnArgs[[2]]
-        s=s+groupsize[i]*mf_tmp
-    }
-}
-ret = s
-return(list( ret= ret,poped.db =poped.db )) 
-}
diff --git a/R/v.R b/R/v.R
index fd6c9409..763f350e 100644
--- a/R/v.R
+++ b/R/v.R
@@ -66,7 +66,7 @@ v <- function(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,d,sigma,docc,poped.db)
         returnArgs <- LinMatrixLH(model_switch,xt_ind,x,a,bpop,b_ind,bocc_ind,size(sigma,1),poped.db) 
         lh <- returnArgs[[1]]
         poped.db <- returnArgs[[2]]
-        interact=lh%*%kron_tmp(d,sigma)%*%t(lh)
+        interact=lh %*% (d %x% sigma) %*% t(lh)
         if (sum(dim(interact))==2){
           ret = ret + interact
         } else {
diff --git a/tests/testthat/examples_fcn_doc/examples_mf.R b/tests/testthat/examples_fcn_doc/examples_mf.R
index 4ce01385..50cd28a4 100644
--- a/tests/testthat/examples_fcn_doc/examples_mf.R
+++ b/tests/testthat/examples_fcn_doc/examples_mf.R
@@ -3,8 +3,9 @@
 ind=1
 
 # no occasion defined in this example, so result is zero
-output <- mf(model_switch=t(poped.db$design$model_switch[ind,,drop=FALSE]),
-   xt_ind=t(poped.db$design$xt[ind,,drop=FALSE]),
+poped.db$settings$iFIMCalculationType = 0
+output <- mf3(model_switch=t(poped.db$design$model_switch[ind,,drop=FALSE]),
+   xt=t(poped.db$design$xt[ind,,drop=FALSE]),
    x=zeros(0,1),
    a=t(poped.db$design$a[ind,,drop=FALSE]),
    bpop=poped.db$parameters$bpop[,2,drop=FALSE],
@@ -12,7 +13,7 @@ output <- mf(model_switch=t(poped.db$design$model_switch[ind,,drop=FALSE]),
    sigma=poped.db$parameters$sigma,
    docc=poped.db$parameters$param.pt.val$docc,
    poped.db)
-
+poped.db$settings$iFIMCalculationType = 1
 # in this simple case the full FIM is just the sum of the individual FIMs
 # and all the individual FIMs are the same
 det(output$ret*32) == det(evaluate.fim(poped.db,fim.calc.type=0))  
diff --git a/tests/testthat/examples_fcn_doc/examples_mf5.R b/tests/testthat/examples_fcn_doc/examples_mf5.R
index 6801334d..667c9861 100644
--- a/tests/testthat/examples_fcn_doc/examples_mf5.R
+++ b/tests/testthat/examples_fcn_doc/examples_mf5.R
@@ -1,9 +1,11 @@
 
 #for the FO approximation
 ind=1
+#for evaluating previous mf5 method
+poped.db$settings$iFIMCalculationType = 4
 
 # no occasion defined in this example, so result is zero
-output <- mf5(model_switch=t(poped.db$design$model_switch[ind,,drop=FALSE]),
+output <- mf3(model_switch=t(poped.db$design$model_switch[ind,,drop=FALSE]),
    xt=t(poped.db$design$xt[ind,,drop=FALSE]),
    x=zeros(0,1),
    a=t(poped.db$design$a[ind,,drop=FALSE]),
diff --git a/tests/testthat/examples_fcn_doc/examples_mf6.R b/tests/testthat/examples_fcn_doc/examples_mf6.R
index a38c14ee..74e680a4 100644
--- a/tests/testthat/examples_fcn_doc/examples_mf6.R
+++ b/tests/testthat/examples_fcn_doc/examples_mf6.R
@@ -3,7 +3,8 @@
 ind=1
 
 # no occasion defined in this example
-output <- mf6(model_switch=t(poped.db$design$model_switch[ind,,drop=FALSE]),
+poped.db$settings$iFIMCalculationType = 5
+output <- mf3(model_switch=t(poped.db$design$model_switch[ind,,drop=FALSE]),
    xt=t(poped.db$design$xt[ind,,drop=FALSE]),
    x=zeros(0,1),
    a=t(poped.db$design$a[ind,,drop=FALSE]),
@@ -12,6 +13,7 @@ output <- mf6(model_switch=t(poped.db$design$model_switch[ind,,drop=FALSE]),
    sigma=poped.db$parameters$sigma,
    docc=poped.db$parameters$param.pt.val$docc,
    poped.db)
+poped.db$settings$iFIMCalculationType = 1
 
 # in this simple case the full FIM is just the sum of the individual FIMs
 # and all the individual FIMs are the same
diff --git a/tests/testthat/examples_fcn_doc/examples_mf8.R b/tests/testthat/examples_fcn_doc/examples_mf8.R
index af329939..158a1011 100644
--- a/tests/testthat/examples_fcn_doc/examples_mf8.R
+++ b/tests/testthat/examples_fcn_doc/examples_mf8.R
@@ -1,9 +1,10 @@
 
 #for the FO approximation
 ind=1
+poped.db$settings$iFIMCalculationType = 7
 
 # no occasion defined in this example, so result is zero
-output <- mf8(model_switch=t(poped.db$design$model_switch[ind,,drop=FALSE]),
+output <- mf3(model_switch=t(poped.db$design$model_switch[ind,,drop=FALSE]),
    xt=t(poped.db$design$xt[ind,,drop=FALSE]),
    x=zeros(0,1),
    a=t(poped.db$design$a[ind,,drop=FALSE]),