backward

clear
p=input("Enter value")
h=input("Enter value")
x=[1981 1991 2091 2191 2291]
y=[363 391 421 467 501]
u=(p-x(1))/h
delay0=y(1)-y(2)
delay1=y(2)-y(3)
delay2=y(3)-y(4)
delay3=y(4)-y(5)
fprintf('delay0=%d\n delay1=%d\n delay2=%d\n delay3=%d\n', delay0, delay1, delay2, delay3)

delsqry0=delay0-delay1
delsqry1=delay1-delay2
delsqry2=delay2-delay3
fprintf('delsqry0=%d\n delsqry1=%d\n delsqry2=%d\n', delsqry0, delsqry1, delsqry2)

delcubey0=delsqry0-delsqry1
delcubey1=delsqry1-delsqry2
fprintf('delcube0=%d\n delcubey1=%d\n ', delcubey0, delcubey1)

delfoury0=delcubey0-delcubey1
fprintf('delfoury0=%d\n  ', delfoury0)

#fx=y(1)+u*delay0+((u*(u-1))/factorial(2))*delsqry0+((u*(u-1)*(u-2))/factorial(3))*delcubey0+((u*(u-1)*(u-2)*(u-3))/factorial(4))*delfoury0
fx=y(1)+u*delay0+((u*(u+1))/factorial(2))*delsqry0+((u*(u+1)*(u+2))/factorial(3))*delcubey0+((u*(u+1)*(u+2)*(u+3))/factorial(4))*delfoury0
fprintf('the populaion for %0.4f', fx)