Perhaps you would benefit from an exploratory tool. Splitting the data into deciles of the x coordinate appears to have been performed in that spirit. With modifications described below, it's a perfectly fine approach.
Many bivariate exploratory methods have been invented. A simple one proposed by John Tukey (EDA, Addison-Wesley 1977) is his "wandering schematic plot." You slice the x-coordinate into bins, erect a vertical boxplot of the corresponding y data at the median of each bin, and connect the key parts of the boxplots (medians, hinges, etc.) into curves (optionally smoothing them). These "wandering traces" provide a picture of the bivariate distribution of the data and allow immediate visual assessment of correlation, linearity of relationship, outliers, and marginal distributions, as well as robust estimation and goodness-of-fit evaluation of any nonlinear regression function.
To this idea Tukey added the thought, consistent with the boxplot idea, that a good way to probe the distribution of data is to start at the middle and work outwards, halving the amount of data as you go. That is, the bins to use need not be cut at equally-spaced quantiles, but instead should reflect the quantiles at the points $2^{-k}$ and $1-2^{-k}$ for $k=1, 2, 3, \ldots$.
To display the varying bin populations we can make each boxplot's width proportional to the amount of data it represents.
The resulting wandering schematic plot would look something like this. Data, as developed from the data summary, are shown as gray dots in the background. Over this the wandering schematic plot has been drawn, with the five traces in color and the boxplots (including any outliers shown) in black and white.
The nature of the near-zero correlation becomes immediately clear: the data twist around. Near their center, ranging from $x=-4$ to $x=4$, they have a strong positive correlation. At extreme values, these data exhibit curvilinear relationships that tend on the whole to be negative. The net correlation coefficient (which happens to be $-0.074$ for these data) is close to zero. However, insisting on interpreting that as "nearly no correlation" or "significant but low correlation" would be the same error spoofed in the old joke about the statistician who was happy with her head in the oven and feet in the icebox because on average the temperature was comfortable. Sometimes a single number just won't do to describe the situation.
Alternative exploratory tools with similar purposes include robust smooths of windowed quantiles of the data and fits of quantile regressions using a range of quantiles. With the ready availability of software to perform these calculations they have perhaps become easier to execute than a wandering schematic trace, but they do not enjoy the same simplicity of construction, ease of interpretation, and broad applicability.
The following R code produced the figure and can be applied to the original data with little or no change. (Ignore the warnings produced by bplt (called by bxp): it complains when it has no outliers to draw.)
#
# Data
#
set.seed(17)
n <- 1449
x <- sort(rnorm(n, 0, 4))
s <- spline(quantile(x, seq(0,1,1/10)), c(0,.03,-.6,.5,-.1,.6,1.2,.7,1.4,.1,.6),
xout=x, method="natural")
#plot(s, type="l")
e <- rnorm(length(x), sd=1)
y <- s$y + e # ($ interferes with MathJax processing on SE)
#
# Calculations
#
q <- 2^(-(2:floor(log(n/10, 2))))
q <- c(rev(q), 1/2, 1-q)
n.bins <- length(q)+1
bins <- cut(x, quantile(x, probs = c(0,q,1)))
x.binmed <- by(x, bins, median)
x.bincount <- by(x, bins, length)
x.bincount.max <- max(x.bincount)
x.delta <- diff(range(x))
cor(x,y)
#
# Plot
#
par(mfrow=c(1,1))
b <- boxplot(y ~ bins, varwidth=TRUE, plot=FALSE)
plot(x,y, pch=19, col="#00000010",
main="Wandering schematic plot", xlab="X", ylab="Y")
for (i in 1:n.bins) {
invisible(bxp(list(stats=b$stats[,i, drop=FALSE],
n=b$n[i],
conf=b$conf[,i, drop=FALSE],
out=b$out[b$group==i],
group=1,
names=b$names[i]), add=TRUE,
boxwex=2*x.delta*x.bincount[i]/x.bincount.max/n.bins,
at=x.binmed[i]))
}
colors <- hsv(seq(2/6, 1, 1/6), 3/4, 5/6)
temp <- sapply(1:5, function(i) lines(spline(x.binmed, b$stats[i,],
method="natural"), col=colors[i], lwd=2))
