First enter the data

xd={Pi/6,Pi/4,Pi/3}
yd=Cos[xd]
ypd=-Sin[xd]

Next, define the Lagrange Polynomials that are used.

L[k_,x_]:=Product[((x-xd[[j]])/(xd[[k]]-xd[[j]])),{j,1,k-1}]*Product[((x-xd[[j]])/(xd[[k]]-xd[[j]])),{j,k+1,3}]

Note: to evaluate the derivative at the data points, use a substitution rule, shown in H.

H[k_,x_]:=(1-2(x-xd[[k]])*D[L[k,xx],xx]/.{xx->xd[[k]]})*L[k,x]^2
HHat[k_,x_]:=(x-xd[[k]])*L[k,x]^2

Now the Hermite Polynomial:

Herman[x_]=Sum[yd[[j]]*H[j,x],{j,1,3}]+Sum[ypd[[j]]*HHat[j,x],{j,1,3}];

To graph the original function and the osculating polynomial:

Plot[{Herman[x],Cos[x]},{x,0,Pi/2},PlotStyle->{Directive[Blue,Dashed,Thick],Red},PlotRange->{-.5,1.5}]

It’s hard to see the difference, so plot the absolute error:

Plot[Abs[Herman[x]-Cos[x]],{x,0,Pi/2},Ticks->{xd,Automatic}]

Notice that in the range of data the error is smallest:

Plot[Abs[Herman[x]-Cos[x]],{x,.5,Pi/3},Ticks->{xd,Automatic}]

Finally, put the polynomial in standard form (with decimal coefficients):

Simplify[Together[Expand[N[Herman[x]]]]]