Best way (program/statistical DSL/library) to curve fit data containing multiple y-values for some x-values My data is as follows:
X1   Y2         Y3         Y4         Y5     Y6
0    1          1          1          1      1  
1    1.33333    -          1.22222    -      1.2
2    -          1.90909    2          -      1.8
3    -          -          2.6        -      2.2
4    2          -          3.51111    2.1    3.2
5    -          -          4          -      4.6
6    2.33333    3.09091    5.28889    2.9    5.6
7    2.53333    2.81818    5.37778    3      5.8
8    2.23333    3.27273    -          3.2    6.4
9    2.56666    4          -          4.1    8.4


I currently have the above data. As you can see I have multiple Y values for some X values. This is because each Y column represents data from a different source. What I want is for a program (like SciDavis, which I am currently using) to curve fit a function to all of this data, not just one X-Y column combination like SciDavis is doing when I (in other words I want to see a nice curve in between the light blue and the red/dark blue graphs, an average curve so the speak).
I have tried numerous search terms to find a way to do this in SciDavis, and to find other programs to do this, but to no avail. I am hoping that someone experienced in this can give me a good suggestion for a program that can do this (or a domain specific programming language/statistical package combined with a library if need be. I can program as well, in fact I am trying to find a function that I can use as a part of a model that I will be using in a simulation that I am coding)
 A: Use any regression analysis program.
Organize your data into two columns, one for x and one for y. You probably also want to include a third column that denotes the column of y that each value came from.
For each value of x you will have several rows, one for each observed y.
Estimate a model using a polynomial, probably in x and x^2.
If you want a regression for specific subsets of the original y columns then restrict your analysis to those cases.
A: It would be best to explain more about what you're doing and what you're hoping to achieve.
Taking it at face value, as mentioned in comments, I'd try a different scale -- specifically, since all the data start from exactly 1 and "spread out" from there -- I'd subtract 1 and take logs from the y-values, and if smoothness is expected, then I'd fit a smooth curve against log(x) (with natural splines at low df, say), and use the y-column indicator in the model. Then estimating the overall fit at the average level, and transforming everything back (exponentiate and add 1). The result would be like a geometric-mean curve (of the data-1).

If you wanted something nearer the arithmetic mean you'd add an upshift of $s^2/2$ before the exponentiation (where the level shift is), which would just tilt it up a little -- but I wouldn't do that in this case; the higher curves are "noisier" on the raw scale and I think the geometric mean is the better choice.
[If you don't seek a smooth curve you just fit a factor in levels of X as well as the column indicator.]
Here's how I did that:
pel=read.table(stdin(),head=TRUE)
    X1   Y2         Y3         Y4         Y5     Y6
    0    1          1          1          1      1  
    1    1.33333    NA         1.22222    NA     1.2
    2    NA         1.90909    2          NA     1.8
    3    NA         NA         2.6        NA     2.2
    4    2          NA         3.51111    2.1    3.2
    5    NA         NA         4          NA     4.6
    6    2.33333    3.09091    5.28889    2.9    5.6
    7    2.53333    2.81818    5.37778    3      5.8
    8    2.23333    3.27273    NA         3.2    6.4
    9    2.56666    4          NA         4.1    8.4

pels=(data.frame(X=rep(0:9,5),stack(pel[,-1])) ) 
tpels=data.frame(lvm1=log(pels$values-1),lX=log(pels$X),ind=pels$ind)[pels$values>1,]

library(splines)
lfit4=lm(lvm1~ns(lX,4)+ind,tpels)
mlev4=sum(lfit4$coefficients[6:9])/5
newp=data.frame(lX=log(seq(.1,9,.1)),ind=rep(tpels$ind[1],90))
ppred=predict(lfit4,newdata=newp)
plot(values~X,pels,col=ind,pch=16)
lines((exp(ppred+mlev4)+1)~exp(newp$lX),col="dimgrey",lwd=2)

