Quick guide to Mathematica

To get help, place the cursor on the function and press F1
To run all commands in one cell, place the cursor there and press Shift + Enter or Enter on the numerical keyboard.
Output of many commands is suppressed by ; at the end of the command. If you want to see the output, delete ;
Argument of functions must be in square brackets [...], braces {...} are for arrays, ranges etc., use double brackets [[...]] to access elements of arrays
expr /. {x -> a} means substitute a for x in expr, expr // function means apply function to expr, i.e. it is equivalent to function[expr].
D[f[x],x] = derivative of f[x] with respect to x.
To define a function of x one can use f[x_] := 1+(x^2), notice the underscore and the colon, or # nad & like in f := 1+((#)^2) &
N[expr] evaluates expr with machine precision, if you use the decimal point in the expression it will be also evaluated with machine precision
Special characters can be inserted by pressing Esc ... Esc, e.g. Esc p Esc gives π
Clear[...] is used to unset any variables which could have been assigned previously.
Some other useful commands: Simplify, Expand, Factor; Integrate, Series, Sum

Clear all symbols from previous evaluations to avoid problems

Clear["Global`*"]

Numerical differentiation

How to find difference approximations to the derivatives?

Direct differentiation of interpolating polynomials. Notice that Mathematica uses the Newton interpolating polynomial, which is, of course, the same polynomial as the Lagrange interpolating polynomial, but it is constructed in a different way

(Subscript[P, n])(x)=f[(Subscript[x, 0])]+(x-(Subscript[x, 0]))f[(Subscript[x, 0]),(Subscript[x, 1])]+⋯+(x-(Subscript[x, 0]))(x-(Subscript[x, 1]))…(x-(Subscript[x, n-1]))f[(Subscript[x, 0]),(Subscript[x, 1]),…,(Subscript[x, n])]

where

f[(Subscript[x, 0])]=f((Subscript[x, 0])),   f[(Subscript[x, 0]),(Subscript[x, 1])]=((f[(Subscript[x, 0])]-f[(Subscript[x, 1])])/((Subscript[x, 0])-(Subscript[x, 1]))),   ⋯,   f[(Subscript[x, 0]),(Subscript[x, 1]),…,(Subscript[x, n])]=((f[(Subscript[x, 0]),(Subscript[x, 1]),…,(Subscript[x, n-1])]-f[(Subscript[x, 1]),(Subscript[x, 2]),…,(Subscript[x, n])])/((Subscript[x, 0])-(Subscript[x, n])))

np=3; (* number of points *)
z=x0; (* point at which we want to find the derivative *)
d=1;  (* order of the derivative *)
grid={x0,x1,x2,x3,x4};
fgrid=Map[f, grid];
points=Table[{grid[[i]],fgrid[[i]]},{i,1,np}];
Print["Used points:", points]
P[x_]=InterpolatingPolynomial[points,x]
myder=D[P[x],{x,d}]/.{x->z};
myder/.Table[grid[[i+1]]->x0+i*h,{i,1,np-1}]//Simplify

Taylor series give formulas and error terms for np-point difference formula, all derivative of order less or equal np - 1  will be determined

np=3; (* number of points *)
grid={x0,x1,x2,x3,x4};
Taylor=Array[1,np];
TaylorEq=Array[1,np];
Do[
Taylor[[i]]=Normal[Series[f[grid[[i]]],{grid[[i]],x,np}]];
TaylorEq[[i]]=f[grid[[i]]]==Taylor[[i]],
{i,1,np}
]
Column[TaylorEq]
unknowns=Flatten[{f[x],Table[D[f[x],{x,i}],{i,1,np-1}]}]
sol=Solve[TaylorEq,unknowns]//Simplify;
z=x0;
resGenGrid=sol/.{x->z}//Simplify
substitutions=Table[grid[[i+1]]->x0+i*h,{i,1,np-1}]
resEquiGrid=(sol/.{x->z}//Simplify)/.substitutions//Simplify

To get coefficients for derivatives from the results above, one can use the built-in function Coefficient

Print["General grid:"]
Do[
Print["Derivative of order ",i,": ",Coefficient[Delete[((D[f[x],{x,i}]/.{x->z})/.resGenGrid),0],Table[f[grid[[i]]],{i,1,np}]]//Simplify],
{i,1,np-1}
]
Print["Equidistant grid:"]
Do[
Print["Derivative of order ",i,": ",Coefficient[Delete[(D[f[x],{x,i}]/.{x->z})/.resEquiGrid,0],Table[f[x0+i*h],{i,0,np-1}]]],
{i,1,np-1}
]

Build-in function for derivatives at all points with a preset order (see help on Numerical Solution of Partial Differential Equations, section Spatial Derivative Approximations for more information about this function)

np=10;
grid = h Range[0,np];
fgrid=Map[f, grid];
coef=NDSolve`FiniteDifferenceDerivative[Derivative[2],grid, fgrid, "DifferenceOrder"->2]

What is the optimal choice of h to obtain the most accurate results?

There is a different optimal h for different difference formulas:
D1(h), D2(h) etc are difference formulas with leading error terms given by e1(h), e2(h) etc.
Comparison of numerical error with functions e1(h), e2(h) etc. can be seen in the figure. 
Round-off error should be less than 2ϵ|f((Subscript[x, 0]))|/h, plotted as a function eroff(h)

n=40;
hg=Table[(2.0^(-i)),{i,1,n}]; (* grid of h*)
ϵ=(2.0^(-53)); (* machine epsilon *)
z=1.0; (* point at which we want to compute the derivative *)
g[x_]=Sin[x]; (* function and exact derivatives needed for error terms *)
exactDer=Table[Abs[D[g[x],{x,i}]/.{x->z}],{i,1,7}];
(* basic forward 2-point difference formula should be O(h) *)
D1[h_]=((-g[z]+g[h+z])/h); 
e1[h_]=exactDer[[2]] h/2;
(* centered 3-point difference formula should be O((h^2)) *)
D2[h_]=((-g[-h+z]+g[h+z])/(2 h));
e2[h_]=exactDer[[3]] (h^2)/6;
(* centered 5-point difference should be O((h^4)) *)
D3[h_]=((g[-2 h+z]-8 g[-h+z]+8 g[h+z]-g[2 h+z])/(12 h));
e3[h_]=exactDer[[5]] (h^4)/30;
(* centered 7-point difference should be O((h^6)) *)
D4[h_]=(1/(60 h))(-g[-3 h+z]+9 g[-2 h+z]-45 g[-h+z]+45 g[h+z]-9 g[2 h+z]+g[3 h+z]);
e4[h_]=exactDer[[7]] (h^6)/140;
(* round-off error *)
eroff[h_]=2ϵ Abs[g[z]]/h;
dataD=Table[Table[{hg[[i]],Abs[exactDer[[1]]-f[hg[[i]]]]},{i,1,n}],{f,{D1,D2,D3,D4}}];
datah=Table[Table[{hg[[i]],f[hg[[i]]]},{i,1,n}],{f,{e1,e2,e3,e4,eroff}}];
ListLogLogPlot[Join[dataD,datah],PlotStyle->{Red,Green,Blue,Brown,Black,Black,Black,Black,Black},Joined->{False,False,False,False,True,True,True,True,True},PlotRange->{(10^(-16)),1}]