I worked out a gnuplot-only way to do this, it involves a bit of work and you'll probably have to fine tune the details for your particular problem.
As an example I generated a data file containing values for the function exp(x) and its Taylor expansions from order zero (T^(0)[exp(x)] = 1) to order 3 (T^(3)[exp(x)] = 1 + x + x**2/2. + x**3/6.). This kind of data is suited to this problem because you will have a high data density around the origin, where all the approximations converge to the exact value, and lower data density away from it. It can be generated like this with gnuplot:
set xrange [0:1]
set table
set output "| grep -v '^$' > data"
plot exp(x), 1, 1+x, 1+x+x**2/2., 1+x+x**2/2.+x**3/6.
unset table ; unset output
Note I'm formatting the output so my data file has no blank lines, otherwise gnuplot treats fields separated by blank lines as different data blocks and this eventually messes up the histograms below. This data looks like this (plot "data"):
Now, I create a 2D histogram with this data. It would be extremely helpful if gnuplot offered this feature, but it doesn't, so the task gets a bit tricky. What I will do is create several 1D histograms. For more info on how to generate the latter, check this.
The first thing is to figure out the width along x and y for your bins, xwidth and ywidth, where the number of data points are counted, that is, we divide the data space into a grid where each element measures xwidth by ywidth and is assigned a number equal to the number of data points contained within. The smaller these elements the better resolution your graph will have, but also the more data points you'll need for it to look good. For my data above, this could be something like
xwidth = 0.02
ywidth = 0.05
Now we declare a function to define our 1D bins (details):
bin(x,width)=width*floor(x/width)+width/2.0
and define the number of bins along each direction. Because the xrange for my data is [0:1] and my yrange is [1:2.8], the number of bins would be 50 and 36, respectively. I could use Nx = xrange / xwidth but that would lead to a float Nx and I want an integer. To be safe I do:
Nx = 50
Ny = 36
It might make more sense to define these values the other way around: calculate xwidth as xrange / Nx, in which case you should not have problems with integer/float.
Now I generate the 1D histograms along y, looping over x values:
set output "| grep -v 'u\\|^$' | sed 's/#/\\n#/g' > data2"
set table
plot for [i=0:(Nx-1)] "./data" using \
(bin($2,ywidth)):( i*xwidth <= $1 && (i+1.)*xwidth > $1 ? 1.0 : 0.0) \
smooth freq
unset table ; unset output
Now data2 contains Nx blocks of data, each of them being a scan along y with Ny data points. The value of these data points is the number of data entries in the original data file. As it is, data2 contains 2D data (y, color), which I need to remap to 3D. The x value is given by the data block position, accessible with the every option in gnuplot. To plot this 3-dimensionally I do:
set output "| grep -v 'u\\|^$' | sed 's/#/\\n#/g' > data3"
set table
splot for [i=0:(Nx-1)] "./data2" every :::i::i using \
((i+0.5)*xwidth):1:2
unset table ; unset output
This data3 can now be plotted as a color map:
plot "./data3" with image
which looks like this:
Had I used higher quality data (i.e. with higher resolution) the graph would look nicer. With 2x resolution along each direction, the same looks like below:
