forward

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(2)-y(1)
delay1=y(3)-y(2)
delay2=y(4)-y(3)
delay3=y(5)-y(4)
fprintf('delay0=%d\n delay1=%d\n delay2=%d\n delay3=%d\n', delay0, delay1, delay2, delay3)

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

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

delfoury0=delcubey1-delcubey0
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)