If you need higher-order derivatives, you won't get good results using equidistant data points. If you can sample your function at arbitrary nodes, I would recommend using Chebyshev points, i.e. $$x_k = \cos\left( \pi \frac{k}{n} \right), \quad k=0\dots n$$ for a polynomial of degree $n$.
You can evaluate the polynomial stably using Barycentric Interpolation. Note that since you're using a polynomial of high degree over the entire interval, you will probably need less data points. Note also that this assumes that your data can be represented by a polynomial, i.e. it is continuous and smooth.
Getting higher-order derivatives from the interpolant is a bit tricky and is ill-conditioned for high derivatives. It can be done, however, using Chebyshev polynomials. In Matlab/Octave (sorry, my Python is not good at all), though, you could do the following:
% We will use the sine function as a test case
f = @(x) sin( 4*pi*x );
% Set the number of points and the interval
N = 40;
a = 0; b = 1;
% Create a Vandermonde-like matrix for the interpolation using the
% three-term recurrence relation for the Chebyshev polynomials.
x = cos( pi*[0:N-1]/(N-1) )';
V = ones( N ); V(:,2) = x;
for k=3:N, V(:,k) = 2*x.*V(:,k-1) - V(:,k-2); end;
% Compute the Chebyshev coefficients of the interpolation. Note that we
% map the points x to the interval [a,b]. Note also that the matrix inverse
% can be either computed explicitly or evaluated using a discrete cosine transform.
c = V \ f( (a+b)/2 + (b-a)/2*x );
% Compute the derivative: this is a bit trickier and relies on the relationship
% between Chebyshev polynomials of the first and second kind.
temp = [ 0 ; 0 ; 2*(N-1:-1:1)'.*c(end:-1:2) ];
cdiff = zeros( N+1 , 1 );
cdiff(1:2:end) = cumsum( temp(1:2:end) );
cdiff(2:2:end) = cumsum( temp(2:2:end) );
cdiff(end) = 0.5*cdiff(end);
cdiff = cdiff(end:-1:3);
% Evaluate the polynomial and its derivative at a set of points and plot them.
xx = linspace(-1,1,200)';
Vxx = ones( length(xx) , N ); Vxx(:,2) = xx;
for k=3:N, Vxx(:,k) = 2*xx.*Vxx(:,k-1) - Vxx(:,k-2); end;
plot( (a+b)/2 + (b-a)/2*xx , [ Vxx*c , Vxx(:,1:N-1)*cdiff ] );
The code to compute the derivative can be re-applied several times to compute higher derivatives.
If you use Matlab, you may be interested in the Chebfun project, which does most of this automatically and from which parts of the code example above were taken. Chebfun can create an interpolant from literally any function, e.g. continuous, discontinuous, with singularities, etc... and then compute its integral, derivatives, use it to solve ODEs, etc...
