There are many formulas which might satisfy your requiremements.
Nested powers
One possibility:
$points = 10000 * pow(0.993575964272119, pow($rank, 3.16332422407427) - 1)
This gives you the following results:
f(1) = 10000
f(2) = 9500
f(4) = 6000
f(9) = 12.065
f(10) = 0.84341
f(200) = 0
So the three values you fixed (1, 2 and 4) are all satisfied, but the result for 200 indicates that this might not be exactly what you're looking for. The curve looks like this:

By the way, I found this using python and mpmath, by fixing the form of the formula and determining the numbers with the many digits numerically:
>>> import mpmath
>>> print(mpmath.findroot((lambda a,b: 10000*a**(2**b - 1) - 9500,
... lambda a,b: 10000*a**(4**b - 1) - 6000),
... (0.995, 2.7)))
[0.993575964272119]
[ 3.16332422407427]
If you decide on a different form of the function, this approach might be adapted.
Exp of a polynomial
A possible different form with the desired properties would be this:
$points = exp(9.14265175282929 + $rank*(0.127179575914116 - $rank*0.0594909567672230))
This does not decrease quite as quickly as the one above:
f( 1) = 10000
f( 2) = 9500
f( 4) = 6000
f( 13) = 2.1002
f( 14) = 0.47852
f(200) = 0
It was obtained by solving this system of equations:
a + b + c = log(10000)
a + 2b + 4c = log( 9500)
a + 4b + 16c = log( 6000)
to obtain the coefficients a through c for the polynomial. One can add another degree to match f(200)=2
as well, but in that case, the last coefficient will become positive, which means that points will start to increase with rank for very large ranks.
If you want to match that f(200)=2
as well, you can do so using
$points = exp(max(8.86291000469285 - $rank*0.0408488141206645,
9.14265175282929 + $rank*(0.127179575914116 - $rank*0.0594909567672230)))
although this will result in a bend in your curve.
To compare these alternatives to the above:
