How to tell how small `h` is the in the finite difference?

I'm doing a project required to compute the precise(up to a certain precession N) numerical derivatives of a function. The usual approach was to use the finite difference types of algorithms, i.e.

(f(x+h)-f(x))/h 

However, there's not many document to mention how small the h should be, or how to check the convergence or the stability of the results, as h=0 had a different meaning in the numerical computation than symbolic computation.

How to tell how small h is in the finite difference for the precession N? Also, is there any numerical differentiation method that's stable under recurrent application, since the aim was to apply the differentiation up to 1000 times?

In IEEE754 double precision the optimal value of h depends on f(x) and the value of f''(x), the second derivative of f at x. If you have an upper bound M of its absolute value, then the optimal h is (Scilab code below)

h_opt = 2*sqrt(%eps*f(x)/M)

For example, if f(x)=cos(x) and try to approximate f'(1), you have M=1 and

--> h_opt = 2*sqrt(%eps*cos(1))

h_opt =

   2.191D-08

If you try values of h between 10^-1 and 10^-16 and compute the difference between the forward approximation and the true value -sin(1)

format e
h=10.^-(1:16)';
[h -sin(1)-(cos(1+h)-cos(1))./h]

ans =

   1.000D-01   2.559D-02
   1.000D-02   2.687D-03
   1.000D-03   2.700D-04
   1.000D-04   2.701D-05
   1.000D-05   2.702D-06
   1.000D-06   2.701D-07
   1.000D-07   2.806D-08
   1.000D-08  -3.025D-09
   1.000D-09   1.302D-07
   1.000D-10   3.522D-07
   1.000D-11   3.522D-07
   1.000D-12   7.807D-05
   1.000D-13  -1.032D-03
   1.000D-14   2.299D-03
   1.000D-15   4.671D-02
   1.000D-16  -8.415D-01

you can see that the best precision is attained for h=10^-8, which is well predicted by the magnitude of h_opt. This video https://youtu.be/yeFJWtSnV7Y?t=2013 (in french) explains the computation of the optimal h.