What's the equation behind this strange curve? I have come across a strange kind of curve, which I can't imagine what is the equation behind. I've tried a variety of approaches, from numerically calculating the derivative and second (and so on) derivative... but that gave me no clue about the very nature of this curve. 

@user21820: the data can be found on pastebin: here
This sequence of points seems to grow with a slope of exactly 1 until the slope decreases to zero quite sharply... I numerically calculated the first three derivatives of this curve:
The first derivative looks like :

interestingly (I think), the first derivative's "turning point" is at x=500, just like the y maxima of the function I a looking for... This is numerically verified by the look of the second derivative:

the third derivative, if of any help? (note that the grey line here should not to be considered)

Well... Hope you can help me !
NB: I can provide the raw data if necessary, but I felt this would render the question unreadable (2000 values...)
edit: @Ron: The plot of $f(x) = 500 * \tanh(x/500)$ in red; and $f(x) = 500 * \sqrt(\tanh((x/500)^2))$ in green:

edit2: @user21280: This is also very close but there seems to be a little disagreement, again around the bending point.

edit3: The story behind this mysterious curve:
In a closed system with two compartments A and B, exchanges of some element E occurs between these two compartment. A is finite while B is not. Flux from A to B (Fab) in one unit of time is a fraction of stock of E in A (Fab = 0.001 * A). Flux from B to A is a fraction of stock of E in B times the fraction of free "space" in A (Fba = 0.1 * B * (Amax-A)/Amax
Considering, at start of simulation, that A = 0 and B = x, for x = 1:2000, what is the value y of B at the equilibrium ? (given Amax = 500)

When B0 gets high enough, A is saturated and B_eq (y) increases as a linear function of B0. Hence the idea of substracting B0 from B_eq. Which gives away my mysterious curve.
edit4: the reproducible 'R' code behind:
Amax <- 500 # max stock of element in solid phase

A0  <- 0    # stock of element in solid phase (A) at time zero

niequ <- 1000   # time allowed for reaching equilibrium between A and B

rec <- c(0)

for(B0 in 1:2000){  # B0 = stock of element in B at time zero
    A  <- A0

    B <- B0

    for(i in 1:niequ){

        Fab <- .001*A
        Fba <- 0.1*B*(1-A/Amax)
        B <- B - Fba + Fab
        A <- A  + Fba - Fab

    }
    rec <- c(rec,mean(B))
}
rec <- rec[-1]
plot(rec~c(1:length(rec)))  # B at equilibrium as a function of B at time zero

rec2 <- (rec-c(1:length(rec)))  # substraction of linear increase
rec2 <- -rec2 # for working with positive values (easier for me)

plot(x=c(1:length(rec)), y=rec2)    # the mysterious curve...

 A: This to me seems like a hyperbolic tangent function:
$$f(x) = a \tanh{b x}$$
A quick inspection reveals that $a$ should be about $500$.  Perhaps the slope at $x=0$ being about $1.0$ would imply that $b=1/500$.  Obviously, these numbers should be taken with a grain of salt but at least represent a starting point for the analysis.  The functional form seems to at least match somewhat the derivative profiles.
EDIT
A comparison to the data showed that $f$ did not bend fast enough.  Another possibility is
$$g(x) = 500 \sqrt{\tanh{\left [\left (\frac{x}{500} \right )^2\right ]}}$$
A quick inspection shows that there is a better fit to the bending, although the behavior of the second derivative is a bit different.
A: I don't know a good reason for the data, but it looks like $x \mapsto \int_0^x \left( \frac{1}{2} - \frac{1}{2}\tanh(\frac{t-500}{100}) \right)\ dt$. It matches all the graphs perfectly as far as I can tell.
[Edit: After Rodolphe finally told the origin of the data...]
It can be seen that it is just a differential equation $\frac{da}{dt} = -pa+q(1-\frac{a}{m})(c-a)$ where $a$ is the amount in $A$, $c$ is the total amount in both, $m$ is the capacity of $A$ and $p,q$ are fixed parameters. The differential equation is separable and just requires integrating the reciprocal of a quadratic function, which is easy. I'm too lazy to do it though.
[Edit: After Rodolphe posted the code, it is clear that I misunderstood the process...]
Michael did it so I don't have to. Great!
A: At equilibrium, 
$$F_{ab}=F_{ba}\\
0.001A = 0.1B(1-A/A_{max})\\
0.001A = 0.1(x-A)(1-A/A_{max})\\
0.01AA_{max}=(A-x)(A-A_{max})\\
A^2-(x+1.01A_{max})A+xA_{max}=0\\
A_{eq} = \frac12\left(x+1.01A_{max}\pm\sqrt{x^2-1.98A_{max}x+1.0201A_{max}^2}\right)\\
A_{eq} = \frac12\left(x+505\pm\sqrt{x^2-990x+255025}\right)\\
\text{(probably the negative square-root)}\\
y=x-A_{eq} = \frac12\left(x-1.01A_{max}\mp\sqrt{x^2-1.98A_{max}x+1.0201A_{max}^2}\right)\\
y=\frac12\left(x-505\mp\sqrt{x^2-990x+255025}\right)\\
\text{(probably the positive square-root)}
$$
A: Here is IMHO a slightly more direct approach.
Indeed, this curve with 2 neat asymptotes invites to consider it as a (branch of) hyperbola having a horizontal asymptote with approx. equation $y-500=0$ and a slant asymptote with approx. equation $y-x=0$. The common equation of all hyperbolas with such asymptotes is:
$$(y-500)(y-x)=k$$
where $k$ is constant (where the LHS is the product of the LHS of the asymptotes).
Forcing the curve to pass by a certain point $A(x_0,y_0)$ clearly fixes the value of the constant. If we fix (approx.) $A(1,0.9900796)$ (according to the given data), we obtain the following equation:
$$(y-500)(y-x)=k_0\tag{1}$$
If we expand (1) as a polynomial in $y$ with parameter $x$, we get a quadratic equation
$$y^2-(x+500)y+(500x-k_0)=0$$
whose solutions are:
$$y=\tfrac12 \left(x+500 \pm \sqrt{(x+500)^2-(500x-k_0)} \right)$$
finding back the solution given by others.
Remark: in fact the $\pm$ sign accounts for the two branches of the hyperbola. Clearly, we need to replace it by a minus sign in order to get the desired branch.
