These are just some quick-and-dirty suggestions. Hopefully somebody will find them helpful!

1. Do you have a symbolic function or a set of points?

• If you have a symbolic function, you may be able to calculate the derivative analytically. (Chances are, you would have done this if it were that easy, and you would not be here looking for alternatives.)
• If you have a symbolic function and cannot calculate the derivative analytically, you can always evaluate the function on a set of points, and use some other method listed on this page to evaluate the derivative.
• In most cases, you have a set of points (xi,fi), and will have to use one of the following methods....

2. Is your grid evenly or unevenly spaced?

• If your grid is evenly spaced, you probably will want to use a finite difference scheme (see either of the Wikipedia articles here or here), unless you are using periodic boundary conditions (see below). Here is a decent introduction to finite difference methods in the context of solving ordinary differential equations on a grid (see especially slides 9-14). These methods are generally computationally efficient, simple to implement, and the error of the method can be simply estimated as the truncation error of the Taylor expansions used to derive it.
• If your grid is unevenly spaced, you can still use a finite difference scheme, but the expressions are more difficult and the accuracy varies very strongly with how uniform your grid is. If your grid is very non-uniform, you will probably need to use large stencil sizes (more neighboring points) to calculate the derivative at a given point. People often construct an interpolating polynomial (often the Lagrange polynomial) and differentiate that polynomial to compute the derivative. See for instance, this StackExchange question. It is often difficult to estimate the error using these methods (although some have attempted to do so: here and here). Fornberg's method is often very useful in these cases....
• Care must be taken at the boundaries of your domain because the stencil often involves points that are outside the domain. Some people introduce "ghost points" or combine boundary conditions with derivatives of different orders to eliminate these "ghost points" and simplify the stencil. Another approach is to use right- or left-sided finite difference methods.
• Here's an excellent "cheat sheet" of finite difference methods, including centered, right- and left-sided schemes of low orders. I keep a printout of this near my workstation because I find it so useful.

3. Is your domain periodic? Can you assume periodic boundary conditions?

• If your domain is periodic, you can compute derivatives to a very high order accuracy using Fourier spectral methods. This technique sacrifices performance somewhat to gain high accuracy. In fact, if you are using N points, your estimate of the derivative is approximately N^th order accurate. For more information, see (for example) this WikiBook.
• Fourier methods often use the Fast Fourier Transform (FFT) algorithm to achieve roughly O(N log(N)) performance, rather than the O(N^2) algorithm that a naively-implemented discrete Fourier transform (DFT) might employ.
• If your function and domain are not periodic, you should not use the Fourier spectral method. If you attempt to use it with a function that is not periodic, you will get large errors and undesirable "ringing" phenomena.
• Computing derivatives of any order requires 1) a transform from grid-space to spectral space (O(N log(N))), 2) multiplication of the Fourier coefficients by their spectral wavenumbers (O(N)), and 2) an inverse transform from spectral space to grid space (again O(N log(N))).
• Care must be taken when multiplying the Fourier coefficients by their spectral wavenumbers. Every implementation of the FFT algorithm seems to have its own ordering of the spectral modes and normalization parameters. See, for instance, the answer to this question on the Math StackExchange, for notes about doing this in MATLAB.

4. What level of accuracy are you looking for? Do you need to compute the derivatives within a given tolerance?

• For many purposes, a 1st or 2nd order finite difference scheme may be sufficient. For higher precision, you can use higher order Taylor expansions, dropping higher-order terms.
• If you need to compute the derivatives within a given tolerance, you may want to look around for a high-order scheme that has the error you need.
• Often, the best way to reduce error is reducing the grid spacing in a finite difference scheme, but this is not always possible.
• Be aware that higher-order finite difference schemes almost always require larger stencil sizes (more neighboring points). This can cause issues at the boundaries. (See the discussion above about ghost points.)

5. Does it matter to you that your derivative is evaluated on the same points as your function is defined?

• MATLAB provides the diff function to compute differences between adjacent array elements. This can be used to calculate approximate derivatives via a first-order forward-differencing (or forward finite difference) scheme, but the estimates are low-order estimates. As described in MATLAB's documentation of diff (link), if you input an array of length N, it will return an array of length N-1. When you estimate derivatives using this method on N points, you will only have estimates of the derivative at N-1 points. (Note that this can be used on uneven grids, if they are sorted in ascending order.)
• In most cases, we want the derivative evaluated at all points, which means we want to use something besides the diff method.

6. Do you need to calculate multiple orders of derivatives?

• One can set up a system of equations in which the grid point function values and the 1st and 2nd order derivatives at these points all depend on each other. This can be found by combining Taylor expansions at neighboring points as usual, but keeping the derivative terms rather than cancelling them out, and linking them together with those of neighboring points. These equations can be solved via linear algebra to give not just the first derivative, but the second as well (or higher orders, if set up properly). I believe these are called combined finite difference schemes, and they are often used in conjunction with compact finite difference schemes, which will be discussed next.
• Compact finite difference schemes (link). In these schemes, one sets up a design matrix and calculates the derivatives at all points simultaneously via a matrix solve. They are called "compact" because they are usually designed to require fewer stencil points than ordinary finite difference schemes of comparable accuracy. Because they involve a matrix equation that links all points together, certain compact finite difference schemes are said to have "spectral-like resolution" (e.g. Lele's 1992 paper--excellent!), meaning that they mimic spectral schemes by depending on all nodal values and, because of this, they maintain accuracy at all length scales. In contrast, typical finite difference methods are only locally accurate (the derivative at point #13, for example, ordinarily doesn't depend on the function value at point #200).
• A current area of research is how best to solve for multiple derivatives in a compact stencil. The results of such research, combined, compact finite difference methods, are powerful and widely applicable, though many researchers tend to tune them for particular needs (performance, accuracy, stability, or a particular field of research such as fluid dynamics).

Ready-to-Go Routines

• As described above, one can use the diff function (link to documentation) to compute rough derivatives between adjacent array elements.

• MATLAB's gradient routine (link to documentation) is a great option for many purposes. It implements a second-order, central difference scheme. It has the advantages of computing derivatives in multiple dimensions and supporting arbitrary grid spacing. (Thanks to @thewaywewalk for pointing out this glaring omission!)

• I used Fornberg's method (see above) to develop a small routine (nderiv_fornberg) to calculate finite differences in one dimension for arbitrary grid spacings. I find it easy to use. It uses sided stencils of 6 points at the boundaries and a centered, 5-point stencil in the interior. It is available at the MATLAB File Exchange here.

Conclusion

The field of numerical differentiation is very diverse. For each method listed above, there are many variants with their own set of advantages and disadvantages. This post is hardly a complete treatment of numerical differentiation.

Every application is different. Hopefully this post gives the interested reader an organized list of considerations and resources for choosing a method that suits their own needs.

This community wiki could be improved with code snippets and examples particular to MATLAB.

diff is simply the difference between two adjacent elements in arr, which is exactly why you lose 1 element for using diff once. For example, 10 elements in an array only have 9 differences.

gradient and del2 are for derivatives. Of course, you can use diff to approximate derivative by dividing the difference by the steps. Usually the step is equally-spaced, but it does not have to be. This answers your question why x is not used in the calculation. I mean, it's okay that your x is not uniform-spaced.

So, why gradient gives us an array with the same length of the original array? It is clearly explained in the manual how the boundary is handled,

The gradient values along the edges of the matrix are calculated with single->sided differences, so that

G(:,1) = A(:,2) - A(:,1);

G(:,N) = A(:,N) - A(:,N-1);

Double-gradient and del2 are not necessarily the same, although they are highly correlated. It's all because how you calculate/approximate the 2nd-oder derivatives. The difference is, the former approximates the 2nd derivative by doing 1st derivative twice and the latter directly approximates the 2nd derivative. Please refer to the help manual, the formula are documented.

Okay, do you really want curvature or the 2nd derivative for each point on arr? They are very different. https://en.wikipedia.org/wiki/Curvature#Precise_definition

You can use diff to get the 2nd derivative from the left. Since diff takes the difference from right to left, e.g. x(2)-x(1), you can first flip x from left to right, then use diff. Some codes like,

x=fliplr(x)
first=x./h
second=diff(first)./h

where h is the space between x. Notice I use ./, which idicates that h can be an array (i.e. non-uniform spaced).

C++ calls it std::adjacent_difference.

int X[8] = {1, 1, 2, 3, 5, 8, 13, 21}; 
int Y[8];

std::adjacent_difference( std::begin(X), std::end(X), std::begin(Y) );

Note that the first element of the destination is a direct copy of the source. The remaining elements are differences.

See it work on Compiler Explorer.

You are looking for the numerical gradient I assume.

t0 = 0;
ts = pi/10;
tf = 2*pi;

t  = t0:ts:tf;
x  = sin(t);
dx = gradient(x)/ts

The purpose of this function is a different one (vector fields), but it offers what diff doesn't: input and output vector of equal length.

gradient calculates the central difference between data points. For an array, matrix, or vector with N values in each row, the ith value is defined by

The gradient at the end points, where i=1 and i=N, is calculated with a single-sided difference between the endpoint value and the next adjacent value within the row. If two or more outputs are specified, gradient also calculates central differences along other dimensions. Unlike the diff function, gradient returns an array with the same number of elements as the input.

This code for the most part does what it intends to do: compute the numerical approximation to the derivative. However, the way you're computing the derivative is slightly incorrect. Basically, for each point in your array, you want to subtract the point to the left with the point to the right and divide by 2*h. If you want to do that, you need to post-multiply the vector y and transform it into a column vector.

Transpose y then take D and multiply this with y:

>> d = D*y.'

d =

     4
     8
    12
    -9

However, I'd like to point out that the first and last entry don't make any sense because of your D matrix:

D =

         0    0.2500         0         0
   -0.2500         0    0.2500         0
         0   -0.2500         0    0.2500
         0         0   -0.2500         0

What you are computing here is the central difference for the numerical derivative, which requires a point to the left and a point to the right of the point you are evaluating in your array. Basically, for the first point, you only have the point to the right where the left point is undefined as you are going out of bounds. Similarly for the last point, you only have a point to the left and the point to the right is also out of bounds as it doesn't exist.

The numerical derivative matrix here assumes that the points going out of bounds are 0. It just turns out that you get the derivative correct at this point which is a fluke (i.e. 2*x at x = 1 is 2). The reason why is because the point at x = 0 does give the output of x.^2 being equal to 0, so the omission of this point in the array acts as if it was in your point array to begin with.

The same thing is required for the second derivative:

>> d2 = D*d

d2 =

    2.0000
    2.0000
   -4.2500
   -3.0000

However, bear in mind that for your previous first derivative result d, the last entry is garbage and so if you computed the derivative of this result, you would be using the incorrect right most term when you look at the third entry of d, and so the last two and also the first two entries of d2 are garbage. It just so happens that the second entry of d2 is correct. This is attributed to x = 0 not actually existing in the previous result, but with the assumption that values out of bounds in your point array are 0, we do get the right result.

You should try this on a longer sequence. For example, try doing it on this:

h = 2;
x = 2:h:20;
y = x.^2;

This defines a sequence going from 2 to 20 in steps of h = 2. This is what we get for x, and d:

x =

    2     4     6     8    10    12    14    16    18    20

d =

     4
     8
    12
    16
    20
    24
    28
    32
    36
   -81

As you can see the first and last entry are meaningless as we don't have a left most point for the first entry and a right most point for the second entry to consider. However, the first point by fluke gives us the correct result as we talked about earlier. In general, the first and last entry should not give you the correct results as we don't have a left or right point to account for the derivative. For the rest of the points, you can see that we are correctly computing the derivative, or 2*x.

Let's take a look at the second derivative, d2:

d2 =

    2.0000
    2.0000
    2.0000
    2.0000
    2.0000
    2.0000
    2.0000
    2.0000
  -28.2500
   -9.0000

It so happens that the first two elements are correct, even though they theoretically shouldn't, but we got it by fluke. However, the last two elements aren't correct.

What you need to take away from this is to remember one thing. Remember that when taking numerical derivatives with the central difference, you have to trim off the first n elements and the last n elements, where n is the order of the derivative you are considering.