diff --git a/DESCRIPTION b/DESCRIPTION
index c4a9bfb..5662cfe 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: pse
 Type: Package
 Title: Parameter Space Exploration with Latin Hypercubes
-Version: 0.4.5
-Date: 2016-03-14
+Version: 0.4.6
+Date: 2016-07-19
 Author: Andre Chalom, Paulo Inacio Knegt Lopez de Prado
 Maintainer: Andre Chalom <andrechalom@gmail.com>
 Depends:
diff --git a/NEWS b/NEWS
index 1d8c8cf..49bd35f 100644
--- a/NEWS
+++ b/NEWS
@@ -1,3 +1,13 @@
+2016-07-19  Andre Chalom
+  * Version is 0.4.6
+  * Improved the plotcv function and fixed plotcv and plotecdf to work with multiple 
+  response variables
+  * Introduced Travis CI
+
+2016-03-14  Andre Chalom
+  * Version is 0.4.5
+  * Added maximum iteration parameter (maxIt) for LHS
+
 2015-10-02  Andre Chalom
   * Version is 0.4.4
   * Added Spearman correlation method for LHScorcorr
diff --git a/R/plotcv.R b/R/plotcv.R
index 526f77f..429053d 100644
--- a/R/plotcv.R
+++ b/R/plotcv.R
@@ -1,22 +1,40 @@
 #' @export
+#' @param quant Maximum quantile to be plotted on the ecdf (used to cut off extreme values in the labels)
 #' @rdname plots
-plotcv <- function(obj, stack = FALSE, index.res = 1, col = index.res, ...) { 
-	if (stack)
-		stop("Unimplemented option: stack!");
-	if (index.res != 1)
-		stop ("Unimplemented option: index.res");
+plotcv <- function(obj, stack = FALSE, index.res = 1:get.noutputs(obj), col = index.res, quant = 0.99, ...) { 
+    opar = par(no.readonly=TRUE)
+    on.exit(par(opar))
 	if (class(obj)!="LHS")
 		stop("The first argument should be of class LHS!");
 	if (get.repetitions(obj)<2)
 		stop("Error in function plotcv: the LHS object must have at least two repetitions!")
-	pointwise <- apply(get.results(obj, FALSE), c(1,2), cv)
-	global <- cv(get.results(obj, TRUE))
-	m <- max(pointwise, 1.05*global)
-	mi <- min(pointwise, global)
-	Ecdf(pointwise, xlim=c(mi, m), xlab="pointwise cv", col=col, ...)
-	abline(v=global, lwd=2, lty=3)
-	if (m > 0.8*max(pointwise)) {pos=2} else {pos=4}
-	text(x=global, y=0.1, label="global cv", pos=pos)
+
+	pointwise <- abs(apply(get.results(obj, FALSE), c(1,2), cv))
+	global <- abs(apply(get.results(obj, TRUE), 2, cv))
+
+    if (stack) {
+        dat = vec(pointwise[,index.res])
+        g = rep(index.res, each=get.N(obj))
+        m <- max(quantile(pointwise[,index.res],quant), 1.05*global)
+   	    mi <- min(pointwise[,index.res], global)
+  	    Ecdf(dat, group = g, xlim=c(mi, m), xlab="pointwise cv", col=col, ...)
+        for (i in index.res)
+    	    abline(v=global[i], lwd=2, lty=3, col=col[i])
+    	if (m > 0.8*max(pointwise)) {pos=2} else {pos=4}
+        text(x=global[1], y=0.1, label="global cv", pos=pos)
+    } else {
+        nl = floor(sqrt(length(index.res)))
+        nc = ceiling(length(index.res)/nl)
+        par(mfrow=c(nl, nc))
+        for (i in index.res) {
+            m <- max(quantile(pointwise[,i],quant), 1.05*global)
+            mi <- min(pointwise[,i], global)
+            Ecdf(pointwise[,i], xlim=c(mi, m), xlab="pointwise cv", ...)
+            abline(v=global[i], lwd=2, lty=3)
+            if (m > 0.8*max(pointwise[,i])) {pos=2} else {pos=4}
+            text(x=global[i], y=0.1, label="global cv", pos=pos)
+        }
+    }
 }
 
 #' Coefficient of Variation
diff --git a/R/plotecdf.R b/R/plotecdf.R
index 7d4b43c..32e43ad 100644
--- a/R/plotecdf.R
+++ b/R/plotecdf.R
@@ -44,12 +44,20 @@
 #' @rdname plots
 #' @import graphics Hmisc
 plotecdf <- function (obj, stack=FALSE, index.res =1:get.noutputs(obj), col=index.res, xlab = NULL, ...) {
+    opar = par(no.readonly=TRUE)
+    on.exit(par(opar))
 	if (is.null (xlab)) xlab = obj$res.names
 	if (stack) {
 		if (length(xlab) > 1) xlab = "obj results"
 		dat <- vec(get.results(obj)[,index.res])
 		g <- rep(index.res, each=dim(obj$res)[1])
 		Ecdf(dat, group=g, col=col, xlab=xlab, ...)
-	} else Ecdf(get.results(obj)[,index.res], xlab=xlab, ...)
+	} else {
+        nl = floor(sqrt(length(index.res)))
+        nc = ceiling(length(index.res)/nl)
+        par(mfrow=c(nl, nc))
+        for (i in index.res)
+            Ecdf(get.results(obj)[,i], xlab=xlab, ...)
+    }
 }
 
diff --git a/man/plots.Rd b/man/plots.Rd
index b14a236..7e42092 100644
--- a/man/plots.Rd
+++ b/man/plots.Rd
@@ -7,7 +7,8 @@
 \alias{plotscatter}
 \title{Uncertainty and Sensitivity Plots.}
 \usage{
-plotcv(obj, stack = FALSE, index.res = 1, col = index.res, ...)
+plotcv(obj, stack = FALSE, index.res = 1:get.noutputs(obj),
+  col = index.res, quant = 0.99, ...)
 
 plotecdf(obj, stack = FALSE, index.res = 1:get.noutputs(obj),
   col = index.res, xlab = NULL, ...)
@@ -32,6 +33,8 @@ all variables identified by different colors.}
 
 \item{col}{An optional vector indicating the colors to be used.}
 
+\item{quant}{Maximum quantile to be plotted on the ecdf (used to cut off extreme values in the labels)}
+
 \item{xlab, ylab}{Labels for the x axis (ecdf) or y axis(prcc). 
 The functions use the name provided in the res.names argument from the LHS function if left blank.}