Well, I guess it's a bit late, but everything looked OK up to where you left off. I think it's better to keep everything in terms of $u$ and then consider that since
$$f(a+hu)=P(u)$$
then on differentiating with respect to $u$,
$$hf^{\prime}(a+hu)=P^{\prime}(u)$$
and
$$h^2f^{\prime\prime}(a+hu)=P^{\prime\prime}(u)$$
Given $x$ we can find $u=(x-a)/h$. If you had a table of Stirling polynomials you could plug them into your expression for $P(u)$ ($f(P)$ in your equation on the bottom of your first page) but lacking that it's perhaps easiest to work your expression from right to left: let $P_0(u)=31.99$. Then $P_1(u)=(u-3)P_0(u)/4-14.23$, and so on until $P(u)=P_4(u)=(u-0)P_3(u)/1+19.96$. Here's some Matlab code:
% Newton.m
x = [50 60 70 80 90];
y = [19.96 36.65 58.81 72.21 94.61];
for k = 2:length(y),
y(k:end) = y(k:end)-y(k-1:end-1);
end
P = y(end);
for k = length(y)-1:-1:1,
P = conv([1 -k+1],P)/k;
P(end) = P(end)+y(k);
end
fprintf('u P(u)\n');
for k = 0:length(x)-1,
fprintf('%d %5.2f\n',k,polyval(P,k));
end
P
Pp = [length(P)-1:-1:0].*P;
Pp = Pp(1:end-1)
Ppp = [length(Pp)-1:-1:0].*Pp;
Ppp = Ppp(1:end-1)
h = x(2)-x(1);
a = x(1);
x0 = 51;
u = (x0-a)/h;
fprintf('f(%.1f) = %f\n',x0,polyval(P,u));
fprintf('f''(%.1f) = %f\n',x0,polyval(Pp,u)/h);
fprintf('f''''(%.1f) = %f\n',x0,polyval(Ppp,u)/h^2);
With output:
u P(u)
0 19.96
1 36.65
2 58.81
3 72.21
4 94.61
P =
1.3329 -10.3692 24.5121 1.2142 19.9600
Pp =
5.3317 -31.1075 49.0242 1.2142
Ppp =
15.9950 -62.2150 49.0242
f(51.0) = 20.316302
f'(51.0) = 0.581084
f''(51.0) = 0.429626
That little table is there to verify that $P(u)$ interpolates properly to produce the correct sequence of ordinates. Then the polynomial $P(u)$ and its derivatives $P^{\prime}(u)$ and $P^{\prime\prime}(u)$ are printed out, leading coefficient first. Then $f(51)$, $f^{\prime}(51)$, and $f^{\prime\prime}(51)$ are displayed.