diff --git a/DESCRIPTION b/DESCRIPTION
index 110780f..0b78940 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: bigDM
 Type: Package
 Title: Scalable Bayesian Disease Mapping Models for High-Dimensional Data
-Version: 0.5.2
-Date: 2023-07-19
+Version: 0.5.3
+Date: 2023-09-08
 Authors@R: 
   c(person(given = "Aritz",
            family = "Adin",
@@ -24,7 +24,7 @@ URL: https://github.com/spatialstatisticsupna/bigDM
 BugReports: https://github.com/spatialstatisticsupna/bigDM/issues
 Depends: R (>= 4.0.0)
 Additional_repositories: https://inla.r-inla-download.org/R/stable
-Imports: crayon, doParallel, fastDummies, foreach, future, future.apply, MASS, Matrix, methods, parallel, RColorBrewer, Rdpack, sf, spatialreg, spdep, stats, utils, rlist
+Imports: crayon, doParallel, fastDummies, foreach, future, future.apply, geos, MASS, Matrix, methods, parallel, RColorBrewer, Rdpack, sf, spatialreg, spdep, stats, utils, rlist
 Suggests: bookdown, INLA (>= 22.12.16), knitr, rmarkdown, testthat (>= 3.0.0), tmap
 RdMacros: Rdpack
 License: GPL-3
diff --git a/NAMESPACE b/NAMESPACE
index 43c1b4c..273f744 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -28,6 +28,7 @@ importFrom(fastDummies,dummy_cols)
 importFrom(foreach,"%dopar%")
 importFrom(foreach,foreach)
 importFrom(future.apply,future_mapply)
+importFrom(geos,geos_prepared_intersects)
 importFrom(grDevices,boxplot.stats)
 importFrom(methods,as)
 importFrom(methods,is)
@@ -35,7 +36,6 @@ importFrom(rlist,list.flatten)
 importFrom(sf,st_as_sf)
 importFrom(sf,st_bbox)
 importFrom(sf,st_centroid)
-importFrom(sf,st_contains)
 importFrom(sf,st_distance)
 importFrom(sf,st_drop_geometry)
 importFrom(sf,st_geometry)
diff --git a/NEWS b/NEWS
index 413d415..40397c3 100644
--- a/NEWS
+++ b/NEWS
@@ -1,3 +1,7 @@
+version 0.5.3
+- bug fixed
+- faster implementation of divide_carto() function
+
 version 0.5.2
 - 'X' argument included to STCAR_INLA() function
 - changes in mergeINLA() function
diff --git a/R/CAR_INLA.R b/R/CAR_INLA.R
index 4e16951..1430b14 100644
--- a/R/CAR_INLA.R
+++ b/R/CAR_INLA.R
@@ -265,8 +265,8 @@ CAR_INLA <- function(carto=NULL, ID.area=NULL, ID.group=NULL, O=NULL, E=NULL, X=
 
                 cat(sprintf("+ Model %d of %d",d,D),"\n")
 
-                Rs <- as(Rs,"TsparseMatrix")
-                Rs.Leroux <- as(Rs.Leroux,"TsparseMatrix")
+                Rs <- as(Rs,"Matrix")
+                Rs.Leroux <- as(Rs.Leroux,"Matrix")
                 # Rs <- inla.as.sparse(Rs)
                 # Rs.Leroux <- inla.as.sparse(Rs.Leroux)
 
@@ -308,8 +308,8 @@ CAR_INLA <- function(carto=NULL, ID.area=NULL, ID.group=NULL, O=NULL, E=NULL, X=
 
                 W <- aux$W
                 n <- nrow(W)
-                Rs <- as(Diagonal(n,colSums(W))-W,"TsparseMatrix")
-                Rs.Leroux <- as(Diagonal(n)-Rs,"TsparseMatrix")
+                Rs <- as(Diagonal(n,colSums(W))-W,"Matrix")
+                Rs.Leroux <- as(Diagonal(n)-Rs,"Matrix")
                 # Rs <- inla.as.sparse(Diagonal(n,colSums(W))-W)
                 # Rs.Leroux <- inla.as.sparse(Diagonal(n)-Rs)
 
@@ -471,5 +471,5 @@ CAR_INLA <- function(carto=NULL, ID.area=NULL, ID.group=NULL, O=NULL, E=NULL, X=
   }
 }
 
-utils::globalVariables(c("inla.as.sparse"))
+# utils::globalVariables(c("inla.as.sparse"))
 utils::globalVariables(c("inla.setOption"))
diff --git a/R/STCAR_INLA.R b/R/STCAR_INLA.R
index 766e263..008242b 100644
--- a/R/STCAR_INLA.R
+++ b/R/STCAR_INLA.R
@@ -215,7 +215,7 @@ STCAR_INLA <- function(carto=NULL, data=NULL, ID.area=NULL, ID.year=NULL, ID.gro
         if(temporal=="rw1") dif <- 1
         if(temporal=="rw2") dif <- 2
         D <- diff(diag(T), differences=dif)
-        Rt <- as(t(D)%*%D, "TsparseMatrix")
+        Rt <- as(t(D)%*%D, "Matrix")
         # Rt <- inla.as.sparse(t(D)%*%D)
 
         ## Define hyperprior distributions ##
@@ -352,8 +352,8 @@ STCAR_INLA <- function(carto=NULL, data=NULL, ID.area=NULL, ID.year=NULL, ID.gro
 
                 cat(sprintf("+ Model %d of %d",d,D),"\n")
 
-                Rs <- as(Rs,"TsparseMatrix")
-                Rs.Leroux <- as(Rs.Leroux,"TsparseMatrix")
+                Rs <- as(Rs,"Matrix")
+                Rs.Leroux <- as(Rs.Leroux,"Matrix")
                 # Rs <- inla.as.sparse(Rs)
                 # Rs.Leroux <- inla.as.sparse(Rs.Leroux)
                 S <- nrow(Rs)
@@ -435,8 +435,8 @@ STCAR_INLA <- function(carto=NULL, data=NULL, ID.area=NULL, ID.year=NULL, ID.gro
                 cat("STEP 2: Fitting global model with INLA (this may take a while...)\n")
 
                 W <- aux$W
-                Rs <- as(Diagonal(S,colSums(W))-W, "TsparseMatrix")
-                Rs.Leroux <- as(Diagonal(S)-Rs, "TsparseMatrix")
+                Rs <- as(Diagonal(S,colSums(W))-W, "Matrix")
+                Rs.Leroux <- as(Diagonal(S)-Rs, "Matrix")
                 # Rs <- inla.as.sparse(Diagonal(S,colSums(W))-W)
                 # Rs.Leroux <- inla.as.sparse(Diagonal(S)-Rs)
 
diff --git a/R/divide_carto.R b/R/divide_carto.R
index 3aa10fe..e3bc844 100644
--- a/R/divide_carto.R
+++ b/R/divide_carto.R
@@ -11,8 +11,9 @@
 #'
 #' @return List of \code{sf} objects with the spatial polygons of each subdomain.
 #'
+#' @importFrom geos geos_prepared_intersects
 #' @importFrom RColorBrewer brewer.pal
-#' @importFrom sf st_as_sf st_contains st_intersects st_set_geometry st_union
+#' @importFrom sf st_as_sf st_set_geometry st_union
 #' @importFrom stats aggregate
 #'
 #' @examples
@@ -45,62 +46,63 @@
 #' @export
 divide_carto <- function(carto, ID.group=NULL, k=0, plot=FALSE){
 
-  ## Transform 'SpatialPolygonsDataFrame' object to 'sf' class
-  carto <- sf::st_as_sf(carto)
+        ## Transform 'SpatialPolygonsDataFrame' object to 'sf' class
+        carto <- sf::st_as_sf(carto)
 
-  ## Construct the grouped 'sf' object by ID.group variable ##
-  carto$temp <- seq(1,nrow(carto))
-  Data <- sf::st_set_geometry(carto, NULL)
+        ## Construct the grouped 'sf' object by ID.group variable ##
+        Data <- sf::st_set_geometry(carto, NULL)
+        carto.group <- stats::aggregate(carto[,"geometry"], list(ID.group=Data[,ID.group]), utils::head)
+        D <- nrow(carto.group)
 
-  carto.group <- stats::aggregate(carto[,"geometry"], list(ID.group=Data[,ID.group]), utils::head)
-  D <- nrow(carto.group)
+        ########################
+        ## Disjoint partition ##
+        ########################
+        group.names <- carto.group$ID.group
+        names(group.names) <- group.names
 
-  ########################
-  ## Disjoint partition ##
-  ########################
-  carto.k0 <- vector("list",D)
-  names(carto.k0) <- carto.group$ID.group
+        carto.k0 <- lapply(group.names, function(x) {
+                aux <- carto[st_set_geometry(carto[,ID.group],NULL)==x,]
+                rownames(aux) <- NULL
+                aux
+        })
 
-  for(i in 1:D){
-    loc <- sf::st_contains(carto.group$geometry[i], carto[,"geometry"])
-    carto.k0[[i]] <- merge(carto, data.frame(loc=unlist(loc)), by.x="temp", by.y="loc")
-    carto.k0[[i]]$temp <- NULL
-  }
+        if(k==0){
+                if(plot) lapply(carto.k0, function(x) plot(x$geometry, main=unique(st_set_geometry(x,NULL)[,ID.group])))
+                return(carto.k0)
+        }
 
-  if(k==0){
-    if(plot) lapply(carto.k0, function(x) plot(x$geometry, main=unique(st_set_geometry(x,NULL)[,ID.group])))
-    return(carto.k0)
-  }
+        ############################################
+        ## Partition including k-order neighbours ##
+        ############################################
+        if(k>0){
 
-  ############################################
-  ## Partition including k-order neighbours ##
-  ############################################
-  if(k>0){
+                carto.k <- vector("list",D)
+                names(carto.k) <- names(carto.k0)
 
-    carto.k <- vector("list",D)
-    names(carto.k) <- names(carto.k0)
+                if(plot) color <- RColorBrewer::brewer.pal(k+2,"Blues")
 
-    if(plot) color <- RColorBrewer::brewer.pal(k+2,"Blues")
+                for(i in 1:D){
+                        aux.carto <- carto.group$geometry[i]
 
-    for(i in 1:D){
-      if(plot) plot(carto.k0[[i]]$geometry, col="lightgrey", main=sort(unique(Data[,ID.group]))[i],
-                    xlim=sf::st_bbox(carto.group$geometry[i])[c(1,3)]*c(0.99,1.01),
-                    ylim=sf::st_bbox(carto.group$geometry[i])[c(2,4)]*c(0.99,1.01))
+                        for(j in 1:k){
+                                # loc <- sf::st_intersects(aux.carto, carto[,"geometry"])
+                                # carto.k[[i]] <- merge(carto, data.frame(loc=unlist(loc)), by.x="temp", by.y="loc")
+                                loc <- geos::geos_prepared_intersects(aux.carto, carto[,"geometry"])
+                                carto.k[[i]] <- carto[loc,]
+                                rownames(carto.k[[i]]) <- NULL
 
-      aux.temp <- carto.k0[[i]]$temp
-      aux.carto <- carto.group$geometry[i]
+                                if(plot){
+                                        plot(carto.k[[i]]$geometry, col=color[j+2], main=sort(unique(Data[,ID.group]))[i],
+                                             xlim=sf::st_bbox(carto.group$geometry[i])[c(1,3)]*c(0.99,1.01),
+                                             ylim=sf::st_bbox(carto.group$geometry[i])[c(2,4)]*c(0.99,1.01))
 
-      for(j in 1:k){
-        loc <- sf::st_intersects(aux.carto, carto[,"geometry"])
-        carto.k[[i]] <- merge(carto, data.frame(loc=unlist(loc)), by.x="temp", by.y="loc")
+                                        plot(carto.k0[[i]]$geometry, col=color[1], add=TRUE)
+                                }
 
-        if(plot) plot(carto.k[[i]]$geometry[!(carto.k[[i]]$temp %in% aux.temp)], col=color[j+2], add=TRUE)
-        aux.temp <- carto.k[[i]]$temp
-        aux.carto <- sf::st_union(carto.k[[i]])
-      }
-      carto.k[[i]]$temp <- NULL
-    }
+                                aux.carto <- sf::st_union(carto.k[[i]])
+                        }
+                }
 
-    return(carto.k)
-  }
+                return(carto.k)
+        }
 }
diff --git a/README.md b/README.md
index a7c3b4b..e0e874d 100644
--- a/README.md
+++ b/README.md
@@ -91,6 +91,10 @@ When using this package, please cite the following papers:
 ```
 news(package="bigDM")
 ```
+__Changes in version 0.5.3__ (2023 Sep 08)
+* small bugs fixed
+* faster implementation of `divide_carto()` function
+
 __Changes in version 0.5.2__ (2023 Jun 14)
 * changes in `mergeINLA()` function
 * 'X' argument included to `STCAR_INLA()` function