diff --git a/R/depth.R b/R/depth.R
new file mode 100644
index 0000000000000000000000000000000000000000..053875ec5c9d224f1efa46a0efcefb92f4256f0c
--- /dev/null
+++ b/R/depth.R
@@ -0,0 +1,7 @@
+#' Depth of a list
+#'
+#' Returns the depth of a list
+#' @title Depth of a list
+#' @param this list
+#' @return integer
+depth <- function(this) ifelse(is.list(this), 1L + max(sapply(this, depth)), 0L)
diff --git a/R/flattenlist.R b/R/flattenlist.R
new file mode 100644
index 0000000000000000000000000000000000000000..cd411c656ded80c172d7850f8c4a74cb4376803b
--- /dev/null
+++ b/R/flattenlist.R
@@ -0,0 +1,24 @@
+#' Flattens list in a single list of data.frames
+#'
+#' Flattens list. Can maybe be made better. It might end up copying data in
+#' memory!? It might change the order of the elements.
+#' @title Flattens list
+#' @param x List to flatten.
+#' @return A flatten list
+flattenlist <- function(x){
+    (n <- depth(x))
+    if(n == 2){
+        # Its fine
+        return(x)
+    }else if(n ==3){
+        unlist(x, recursive=FALSE)
+    }else{
+        morelists <- sapply(x, function(xprime) class(xprime)[1]=="list")
+        out <- c(x[!morelists], unlist(x[morelists], recursive=FALSE))
+        if(sum(morelists)){ 
+            Recall(out)
+        }else{
+            return(out)
+        }
+    }
+}
diff --git a/R/forecastmodel.R b/R/forecastmodel.R
index b1fc9d00b5b454e90e96e40023e092fa698f2446..37ae9775adeb94908b7f6736db10335ca2f3d696 100644
--- a/R/forecastmodel.R
+++ b/R/forecastmodel.R
@@ -212,11 +212,11 @@ forecastmodel <- R6::R6Class("forecastmodel", public = list(
             if(class(L)[1]=="data.frame"){ return(list(L)) }
             if(class(L)[1]!="list"){ stop(pst("The value returned from evaluating: ",input$expr,", was not a matrix, data.frame or a list of them."))}
             if(class(L[[1]])[1]=="matrix"){ return(lapply(L, function(mat){ return(as.data.frame(mat)) })) }
-            return(L)
+            return(flattenlist(L))
         })
-        # Put together in one data.list
-        L <- structure(do.call(c, L), class="data.list")
-        #
+        # Make it a data.list with no subsubelements (it's maybe not a data.list, since it miss "t", however to take subsets etc., it must be a data.list)
+        L <- flattenlist(L)
+        class(L) <- "data.list"
         return(L)
     },
     #----------------------------------------------------------------
diff --git a/R/operator_multiply.R b/R/operator_multiply.R
index 2ce598f0d54a20521110444b2659c834264555a9..fc3e0ca4c02c31b3e89f4ffa14afd34950db9ce6 100644
--- a/R/operator_multiply.R
+++ b/R/operator_multiply.R
@@ -49,35 +49,39 @@
 #' @export
 
 "%**%" <- function(x, y) {
-    if( is.null(dim(y)) ){
-        ## y is not matrix like
-        lapply(x, function(xx) {
-            xx * y
-        })
+    # If any of them is a list: do recursive calls 
+    if( class(x)[1] == "list" ){
+        return(flattenlist(lapply(x, "%**%", y=y)))
+    }else if(class(y)[1] == "list"){
+        return(flattenlist(lapply(y, "%**%", y=x)))
+    }
+
+    # Do the multiplication
+    # If either is just a vector
+    if(is.null(dim(x)) | is.null(dim(y))){
+        return(x * y)
     }else{
-        ## y is matrix like
-        lapply(x, function(xx) {
-            ## Check if different horizon k columns
-            colmatch <- TRUE
-            if (ncol(xx) != ncol(y)) {
-                colmatch <- FALSE
-            }else if(any(nams(xx) != nams(y))){
-                colmatch <- FALSE
-            }
-            if(!colmatch){
-                ## Not same columns, take only the k in both
-                nms <- nams(xx)[nams(xx) %in% nams(y)]
-                xx <- xx[, nms]
-                y <- y[, nms]
-            }
-            ## Now multiply
-            val <- xx * y
-            ## Must be data.frame
-            if( is.null(dim(val)) ){
-                val <- data.frame(val)
-                nams(val) <- nms
-            }
-            return(val)
-        })
+        # Both are matrices
+        # Check if they have different horizon k columns
+        colmatch <- TRUE
+        if (ncol(x) != ncol(y)) {
+            colmatch <- FALSE
+        }else if(any(nams(x) != nams(y))){
+            colmatch <- FALSE
+        }
+        if(!colmatch){
+            # Not same columns, take only the k in both
+            nms <- nams(x)[nams(x) %in% nams(y)]
+            x <- x[, nms]
+            y <- y[, nms]
+        }
+        # Now multiply
+        val <- x * y
+        # Must be data.frame
+        if( is.null(dim(val)) ){
+            val <- data.frame(val)
+            nams(val) <- nms
+        }
+        return(val)
     }
 }