21 files changed, 607 insertions, 0 deletions
diff --git a/1673/CH3/EX3.10/3_10.sce b/1673/CH3/EX3.10/3_10.sce
new file mode 100755
index 000000000..39678ff78
--- /dev/null
+++ b/1673/CH3/EX3.10/3_10.sce
@@ -0,0 +1,40 @@
+//practical interpolation
+//example 3.10
+//page 97
+clc;clear;close;
+x=[0.61 0.62 0.63 0.64 0.65 0.66 0.67];
+y=[1.840431 1.858928 1.877610 1.896481 1.915541 1.934792 1.954237];
+h=0.01//interval between values of x
+c=1;
+for i=1:6
+    d1(c)=y(i+1)-y(i);
+    c=c+1;
+end
+c=1;
+for i=1:5
+    d2(c)=d1(i+1)-d1(i);
+    c=c+1
+end
+c=1;
+for i=1:4
+    d3(c)=d2(i+1)-d2(i);
+    c=c+1;
+end
+c=1;
+for i=1:3
+    d4(c)=d3(i+1)-d3(i);
+    c=c+1;
+end
+d=[d1(1) d2(1) d3(1) d4(1)];
+x0=0.644;
+p=(x0-x(4))/h;
+y_x=y(4);
+y_x=y_x+p*(d1(3)+d1(4))/2+p^2*(d2(2))/2;//stirling formula
+printf(' the value at %f by stirling formula is : %f\n\n',x0,y_x);
+y_x=y(4);
+y_x=y_x+p*d1(4)+p*(p-1)*(d2(3)+d2(4))/2;
+printf(' the value at %f by bessels formula is : %f\n\n',x0,y_x);
+y_x=y(4);
+q=1-p;
+y_x=q*y(4)+q*(q^2-1)*d2(3)/2+p*y(5)+p*(q^2-1)*d2(4)/2;
+printf(' the value at %f by everrets formula is : %f\n\n',x0,y_x);
diff --git a/1673/CH3/EX3.11/3_11.sce b/1673/CH3/EX3.11/3_11.sce
new file mode 100755
index 000000000..bab40737e
--- /dev/null
+++ b/1673/CH3/EX3.11/3_11.sce
@@ -0,0 +1,37 @@
+//practical interpolation
+//example 3.11
+//page 99
+clc;clear;close;
+x=[0.61 0.62 0.63 0.64 0.65 0.66 0.67];
+y=[1.840431 1.858928 1.877610 1.896481 1.915541 1.934792 1.954237];
+h=0.01//interval between values of x
+c=1;
+for i=1:6
+    d1(c)=y(i+1)-y(i);
+    c=c+1;
+end
+c=1;
+for i=1:5
+    d2(c)=d1(i+1)-d1(i);
+    c=c+1
+end
+c=1;
+for i=1:4
+    d3(c)=d2(i+1)-d2(i);
+    c=c+1;
+end
+c=1;
+for i=1:3
+    d4(c)=d3(i+1)-d3(i);
+    c=c+1;
+end
+d=[d1(1) d2(1) d3(1) d4(1)];
+x0=0.638;
+p=(x0-x(4))/h;
+y_x=y(4);
+y_x=y_x+p*(d1(3)+d1(4))/2+p^2*(d2(2))/2;//stirling formula
+printf(' the value at %f by stirling formula is : %f\n\n',x0,y_x);
+y_x=y(3);
+p=(x0-x(3))/h;
+y_x=y_x+p*d1(3)+p*(p-1)*(d2(2)/2);
+printf(' the value at %f by bessels formula is : %f\n\n',x0,y_x);
diff --git a/1673/CH3/EX3.12/3_12.sce b/1673/CH3/EX3.12/3_12.sce
new file mode 100755
index 000000000..c0273f306
--- /dev/null
+++ b/1673/CH3/EX3.12/3_12.sce
@@ -0,0 +1,36 @@
+//practical interpolation
+//example 3.12
+//page 99
+clc;clear;close;
+x=[1.72 1.73 1.74 1.75 1.76 1.77 1.78];
+y=[0.1790661479 0.1772844100 0.1755204006 0.1737739435 0.1720448638 0.1703329888 0.1686381473];
+h=0.01//interval between values of x
+c=1;
+for i=1:6
+    d1(c)=y(i+1)-y(i);
+    c=c+1;
+end
+c=1;
+for i=1:5
+    d2(c)=d1(i+1)-d1(i);
+    c=c+1
+end
+c=1;
+for i=1:4
+    d3(c)=d2(i+1)-d2(i);
+    c=c+1;
+end
+c=1;
+for i=1:3
+    d4(c)=d3(i+1)-d3(i);
+    c=c+1;
+end
+x0=1.7475;
+y_x=y(3);
+p=(x0-x(3))/h;
+y_x=y_x+p*d1(3)+p*(p-1)*((d2(2)+d2(3))/2)/2;
+printf(' the value at %f by bessels formula is : %0.10f\n\n',x0,y_x);
+y_x=y(4);
+q=1-p;
+y_x=q*y(3)+q*(q^2-1)*d2(2)/6+p*y(4)+p*(p^2-1)*d2(2)/6;
+printf(' the value at %f by everrets formula is : %0.10f\n\n',x0,y_x);
+\ No newline at end of file
diff --git a/1673/CH3/EX3.13/3_13.sce b/1673/CH3/EX3.13/3_13.sce
new file mode 100755
index 000000000..c853a5339
--- /dev/null
+++ b/1673/CH3/EX3.13/3_13.sce
@@ -0,0 +1,20 @@
+//example 3.13
+//lagrange's interpolation formula
+//page 104
+clc;clear;close;
+x=[300 304 305 307];
+y=[2.4771 2.4829 2.4843 2.4871];
+x0=301;
+log_301=0;
+poly(0,'x');
+for i=1:4
+    p=y(i);
+    for j=1:4
+        if i~=j then
+            p=p*((x0-x(j) )/( x(i)-x(j)))
+        end
+    end
+    log_301=log_301+p;
+    end
+disp(log_301,'log_301=');
+ 
+\ No newline at end of file
diff --git a/1673/CH3/EX3.14/3_14.sce b/1673/CH3/EX3.14/3_14.sce
new file mode 100755
index 000000000..abc816f29
--- /dev/null
+++ b/1673/CH3/EX3.14/3_14.sce
@@ -0,0 +1,20 @@
+//example 3.14
+//lagrange's interpolation formula
+//page 105
+clc;clear;close;
+y=[4 12 19];
+x=[1 3 4];
+y_x=7;
+Y_X=0;
+poly(0,'y');
+for i=1:3
+    p=x(i);
+    for j=1:3
+        if i~=j then
+            p=p*((y_x-y(j) )/( y(i)-y(j)))
+        end
+    end
+    Y_X=Y_X+p;
+    end
+disp(Y_X,'Y_X=');
+ 
+\ No newline at end of file
diff --git a/1673/CH3/EX3.15/3_15.sce b/1673/CH3/EX3.15/3_15.sce
new file mode 100755
index 000000000..bb61d663c
--- /dev/null
+++ b/1673/CH3/EX3.15/3_15.sce
@@ -0,0 +1,18 @@
+//example 3.15
+//lagrange's interpolation formula
+//page 105
+clc;clear;close;
+x=[2 2.5 3.0];
+y=[0.69315 0.91629 1.09861];
+deff('y=l0(x)','y=(x-2.5)*(x-3.0)/(-0.5)*(-1.0)')
+x=poly(0,'x');
+disp(l0(x),'l0(x)=');
+deff('y=l1(x)','y=((x-2.0)*(x-3.0))/((0.5)*(-0.5))')
+x=poly(0,'x');
+disp(l1(x),'l1(x)=');
+deff('y=l2(x)','y=((x-2.0)*(x-2.5))/((1.0)*(0.5))')
+x=poly(0,'x');
+disp(l2(x),'l2(x)=');
+f_x=l0(2.7)*y(1)+l1(2.7)*y(2)+l2(2.7)*y(3);
+printf(' the calculated value is %f:',f_x);
+printf('\n\n the error occured in the value is %0.9f',abs(f_x-log(2.7)))
+\ No newline at end of file
diff --git a/1673/CH3/EX3.16/3_16.sce b/1673/CH3/EX3.16/3_16.sce
new file mode 100755
index 000000000..65e5222c3
--- /dev/null
+++ b/1673/CH3/EX3.16/3_16.sce
@@ -0,0 +1,20 @@
+//example 3.16
+//lagrange's interpolation formula
+//page 104
+clc;clear;close;
+x=[0 %pi/4 %pi/2];
+y=[0 0.70711 1.0];
+x0=%pi/6;
+sin_x0=0;
+poly(0,'x');
+for i=1:3
+    p=y(i);
+    for j=1:3
+        if j~=i then
+            p=p*((x0-x(j) )/( x(i)-x(j)))
+        end
+    end
+    sin_x0=sin_x0+p;
+    end
+disp(sin_x0,'sin_x0=');
+ 
+\ No newline at end of file
diff --git a/1673/CH3/EX3.17/3_17.sce b/1673/CH3/EX3.17/3_17.sce
new file mode 100755
index 000000000..9c87d06e0
--- /dev/null
+++ b/1673/CH3/EX3.17/3_17.sce
@@ -0,0 +1,21 @@
+//lagrange's interpolation
+//example 3.17
+//page 106
+clc;clear;close;
+x=[0  3 4];
+y=[-12 12 24];
+//1 appears to be one the roots the polynomial
+for i=1:3
+    r_x(i)=y(i)/(x(i)-1);
+end
+deff('y=l0(x)','y=((x-3)*(x-4))/((-3)*(-4))')
+x=poly(0,'x');
+disp(l0(x),'l0(x)=');
+deff('y=l1(x)','y=((x-0)*(x-4))/((3)*(-1))')
+x=poly(0,'x');
+disp(l1(x),'l1(x)=');
+deff('y=l2(x)','y=((x-0)*(x-3))/((4)*(1))')
+x=poly(0,'x');
+disp(l2(x),'l2(x)=');
+disp(l0(x)*r_x(1)+l1(x)*r_x(2)+l2(x)*r_x(3),'f_(x)=');
+disp((x-1)*(l0(x)*r_x(1)+l1(x)*r_x(2)+l2(x)*r_x(3))','the required polynimial is :')
+\ No newline at end of file
diff --git a/1673/CH3/EX3.18/3_18.sce b/1673/CH3/EX3.18/3_18.sce
new file mode 100755
index 000000000..4c185f6a3
--- /dev/null
+++ b/1673/CH3/EX3.18/3_18.sce
@@ -0,0 +1,23 @@
+//error in lagrange's interpolation formula
+//example 3.18
+//page 107
+clc;clear;close;
+x=[2 2.5 3.0];
+y=[0.69315 0.91629 1.09861];
+deff('y=l0(x)','y=(x-2.5)*(x-3.0)/(-0.5)*(-1.0)')
+x=poly(0,'x');
+disp(l0(x),'l0(x)=');
+deff('y=l1(x)','y=((x-2.0)*(x-3.0))/((0.5)*(-0.5))')
+x=poly(0,'x');
+disp(l1(x),'l1(x)=');
+deff('y=l2(x)','y=((x-2.0)*(x-2.5))/((1.0)*(0.5))')
+x=poly(0,'x');
+disp(l2(x),'l2(x)=');
+f_x=l0(2.7)*y(1)+l1(2.7)*y(2)+l2(2.7)*y(3);
+printf(' the calculated value is %f:',f_x);
+err=abs(f_x-log(2.7));
+deff('y=R_n(x)','y=(((x-2)*(x-2.5)*(x-3))/6)');
+est_err=abs(R_n(2.7)*(2/8))
+if est_err>err then
+    printf('\n\n the error agrees with the actual error')
+end
+\ No newline at end of file
diff --git a/1673/CH3/EX3.19/3_19.sce b/1673/CH3/EX3.19/3_19.sce
new file mode 100755
index 000000000..bc7729d63
--- /dev/null
+++ b/1673/CH3/EX3.19/3_19.sce
@@ -0,0 +1,18 @@
+//error in lagrenge's interpolation
+//example 3.19
+//page 107
+clc;clear;close;
+x=[0 %pi/4 %pi/2];
+y=[0 0.70711 1.0];
+deff('y=l0(x)','y=((x-0)*(x-%pi/2))/((%pi/4)*(-%pi/4))')
+x=poly(0,'x');
+disp(l0(x),'l0(x)=');
+deff('y=l1(x)','y=((x-0)*(x-%pi/4))/((%pi/2)*(%pi/4))')
+x=poly(0,'x');
+disp(l1(x),'l1(x)=');
+f_x=l0(%pi/6)*y(2)+l1(%pi/6)*y(3);
+err=abs(f_x-sin(%pi/6));
+deff('y=f(x)','y=((x-0)*(x-%pi/4)*(x-%pi/2))/6');
+if abs(f(%pi/6))>err then
+    printf('\n\n the error agrees with the actual error')
+end
+\ No newline at end of file
diff --git a/1673/CH3/EX3.21/3_21.sce b/1673/CH3/EX3.21/3_21.sce
new file mode 100755
index 000000000..5c73ddf2f
--- /dev/null
+++ b/1673/CH3/EX3.21/3_21.sce
@@ -0,0 +1,32 @@
+//hermite's interpolation formula
+//exammple 3.21
+//page 110
+clc;clear;close;
+x=[2.0 2.5 3.0]
+y=[0.69315 0.91629 1.09861]
+deff('y=f(x)','y=log(x)')
+h=0.0001;
+for i=1:3
+    y1(i)=(f(x(i)+h)-f(x(i)))/h;
+end
+deff('y=l0(x)','y=(x-2.5)*(x-3.0)/(-0.5)*(-1.0)')
+a=poly(0,'x');
+disp(l0(a),'l0(x)=');
+deff('y=l1(x)','y=((x-2.0)*(x-3.0))/((0.5)*(-0.5))')
+a=poly(0,'x');
+disp(l1(a),'l1(x)=');
+deff('y=l2(x)','y=((x-2.0)*(x-2.5))/((1.0)*(0.5))')
+a=poly(0,'x');
+disp(l2(a),'l2(x)=');
+dl0=(l0(x(1)+h)-l0(x(1)))/h;
+dl1=(l1(x(2)+h)-l1(x(2)))/h;
+dl2=(l2(x(3)+h)-l2(x(3)))/h;
+x0=2.7;
+u0=[1-2*(x0-x(1))*dl0]*(l0(x0))^2;
+u1=[1-2*(x0-x(2))*dl1]*(l1(x0))^2;
+u2=[1-2*(x0-x(3))*dl2]*(l2(x0))^2;
+v0=(x0-x(1))*l0(x0)^2;
+v1=(x0-x(2))*l1(x0)^2;
+v2=(x0-x(3))*l2(x0)^2;
+H=u0*y(1)+u1*y(2)+u2*y(3)+v0*y1(1)+v1*y1(2)+v2*y1(3);
+printf(' the approximate value of ln(%0.2f) is %0.6f:',x0,H);
diff --git a/1673/CH3/EX3.22/3_22.sce b/1673/CH3/EX3.22/3_22.sce
new file mode 100755
index 000000000..e68596aec
--- /dev/null
+++ b/1673/CH3/EX3.22/3_22.sce
@@ -0,0 +1,15 @@
+//newton's general interpolation formula
+//example 3.22
+//page 114
+clc;clear;close;
+x=[300 304 305 307];
+y=[2.4771 2.4829 2.4843 2.4871];
+for i=1:3
+    d1(i)=(y(i+1)-y(i))/(x(i+1)-x(i));
+end
+for i=1:2
+    d2(i)=(d1(i+1)-d1(i))/(x(i+2)-x(i));
+end
+x0=301;
+log301=y(1)+(x0-x(1))*d1(1)+(x0-x(2))*d2(1);
+printf(' valure of log(%d) is :%0.4f',x0,log301);
diff --git a/1673/CH3/EX3.23/3_23.sce b/1673/CH3/EX3.23/3_23.sce
new file mode 100755
index 000000000..522722787
--- /dev/null
+++ b/1673/CH3/EX3.23/3_23.sce
@@ -0,0 +1,21 @@
+//example 3.22
+//newton's divided formula
+//page 114
+clc;clear;close
+x=[-1 0 3 6 7];
+y=[3 -6 39 822 1611];
+for i=1:4
+    d1(i)=(y(i+1)-y(i))/(x(i+1)-x(i));
+end
+for i=1:3
+    d2(i)=(d1(i+1)-d1(i))/(x(i+2)-x(i));
+end
+for i=1:2
+    d3(i)=(d2(i+1)-d2(i))/(x(i+3)-x(i));
+end
+for i=1:1
+    d4(i)=(d3(i+1)-d3(i))/(x(i+4)-x(i));
+end
+X=poly(0,'X')
+f_x=y(1)+(X-x(1))*(d1(1))+(X-x(2))*(X-x(1))*d2(1)+(X-x(1))*(X-x(2))*(X-x(3))*d3(1)+(X-x(1))*(X-x(2))*(X-x(3))*(X-x(4))*d4(1)
+disp(f_x,'the polynomial equation is =')
+\ No newline at end of file
diff --git a/1673/CH3/EX3.24/3_24.sce b/1673/CH3/EX3.24/3_24.sce
new file mode 100755
index 000000000..8373fd68c
--- /dev/null
+++ b/1673/CH3/EX3.24/3_24.sce
@@ -0,0 +1,20 @@
+//interpolation by iteration
+//example 3.24
+//page 116
+clc;clear;close;
+x=[300 304 305 307];
+y=[2.4771 2.4829 2.4843 2.4871];
+x0=301;
+for i=1:3
+    d=determ([y(i),(x(i)-x0);y(i+1),(x(i+1)-x0)])
+    d1(i)=d/(x(i+1)-x(i));
+end
+for i=1:2
+    d=determ([d1(i),(x(i+1)-x0);d1(i+1),(x(i+2)-x0)])
+    d2(i)=d/(x(i+2)-x(i+1));
+end
+for i=1:1
+    d=determ([d2(i),(x(i+2)-x0);d2(i+1),(x(i+3)-x0)])
+    d3(i)=d/(x(i+3)-x(i+2));
+end
+printf(' the value of log(%d) is : %f',x0,d3(1))
diff --git a/1673/CH3/EX3.25/3_25.sce b/1673/CH3/EX3.25/3_25.sce
new file mode 100755
index 000000000..eb79715f6
--- /dev/null
+++ b/1673/CH3/EX3.25/3_25.sce
@@ -0,0 +1,25 @@
+//inverse intrpolation
+//example 3.25
+//page  118
+clc;clear;close;
+x=[2 3 4 5];
+y=[8 27 64 125];
+for i=1:3
+    d1(i)=y(i+1)-y(i);
+end
+for i=1:2
+    d2(i)=d1(i+1)-d1(i);
+end
+for i=1:1
+    d3(i)=d2(i+1)-d2(i);
+end
+yu=10;//square rooot of 10
+y0=y(1);
+d=[d1(1) d2(1) d3(1)];
+u1=(yu-y0)/d1(1);
+u2=((yu-y0-u1*(u1-1)*d2(1)/2)/d1(1));
+u3=(yu-y0-u2*(u2-1)*d2(1)/2-u2*(u2-1)*(u2-2)*d3(1)/6)/d1(1);
+u4=(yu-y0-u3*(u3-1)*d2(1)/2-u3*(u3-1)*(u3-2)*d3(1)/6)/d1(1);
+u5=(yu-y0-u4*(u4-1)*d2(1)/2-u4*(u4-1)*(u4-2)*d3(1)/6)/d1(1);
+printf(' %f \n %f \n %f \n %f \n %f \n ',u1,u2,u3,u4,u5);
+printf(' the approximate square root of %d is: %0.3f',yu,x(1)+u5)
+\ No newline at end of file
diff --git a/1673/CH3/EX3.26/3_26.sce b/1673/CH3/EX3.26/3_26.sce
new file mode 100755
index 000000000..6fd19b207
--- /dev/null
+++ b/1673/CH3/EX3.26/3_26.sce
@@ -0,0 +1,26 @@
+//double interpolation 
+//example 3.26
+//page 119
+clc;clear;close;
+y=[0 1 2 3 4];
+x=[0 1 4 9 16;2 3 6 11 18;6 7 10 15 22;12 13 16 21 28;18 19 22 27 34];
+printf(' y\t\n');
+for i=1:5
+    printf('\n%d',y(i));
+end
+printf('\n\n-----------------------------------------------x--------------------------------------------\n');
+printf('0\t     1\t    2\t    3\t   4\t\n');
+printf('--------------------------------------------------------------------------------------------\n');
+for i=1:5
+    for j=1:5
+printf('%d\t',x(i,j));
+end
+printf('\n');
+end
+//for x=2.5;
+for i=1:5
+    z(i)=(x(i,3)+x(i,4))/2;
+end
+//y=1.5;
+Z=(z(2)+z(3))/2;
+printf(' the interpolated value when x=2.5 and y=1.5 is : %f',Z);
+\ No newline at end of file
diff --git a/1673/CH3/EX3.4/3_4.sce b/1673/CH3/EX3.4/3_4.sce
new file mode 100755
index 000000000..e251cfb51
--- /dev/null
+++ b/1673/CH3/EX3.4/3_4.sce
@@ -0,0 +1,37 @@
+//example 3.4
+//interpolation
+//page 86
+clc;clear;close;
+x=[1 3 5 7];
+y=[24 120 336 720];
+h=2//interval between values of x
+c=1;
+for i=1:3
+    d1(c)=y(i+1)-y(i);
+    c=c+1;
+end
+c=1;
+for i=1:2
+    d2(c)=d1(i+1)-d1(i);
+    c=c+1
+end
+c=1;
+for i=1:1
+    d3(c)=d2(i+1)-d2(i);
+    c=c+1;
+end
+
+d=[d1(1) d2(1) d3(1)];
+x0=8;//value at 8;
+pp=1;
+y_x=y(1);
+p=(x0-1)/2;
+for i=1:3
+    pp=1;
+    for j=1:i
+    pp=pp*(p-(j-1))    
+    end
+y_x=y_x+(pp*d(i))/factorial(i);
+end
+printf('value of function at %f is :%f',x0,y_x);
+
diff --git a/1673/CH3/EX3.6/3_6.sce b/1673/CH3/EX3.6/3_6.sce
new file mode 100755
index 000000000..674e12ed9
--- /dev/null
+++ b/1673/CH3/EX3.6/3_6.sce
@@ -0,0 +1,46 @@
+//example 3.6
+//interpolation
+//page 87
+clc;clear;close;
+x=[15 20 25 30 35 40];
+y=[0.2588190 0.3420201 0.4226183 0.5 0.5735764 0.6427876];
+h=5//interval between values of x
+c=1;
+for i=1:5
+    d1(c)=y(i+1)-y(i);
+    c=c+1;
+end
+c=1;
+for i=1:4
+    d2(c)=d1(i+1)-d1(i);
+    c=c+1
+end
+c=1;
+for i=1:3
+    d3(c)=d2(i+1)-d2(i);
+    c=c+1;
+end
+c=1;
+for i=1:2
+    d4(c)=d3(i+1)-d3(i);
+    c=c+1;
+end
+c=1;
+for i=1:1
+    d5(c)=d4(i+1)-d4(i);
+    c=c+1;
+end
+c=1;
+d=[d1(5) d2(4) d3(3) d4(2) d5(1)];
+x0=38;//value at 38 degree
+pp=1;
+y_x=y(6);
+p=(x0-x(6))/h;
+for i=1:5
+    pp=1;
+    for j=1:i
+    pp=pp*(p+(j-1))   
+end
+y_x=y_x+((pp*d(i))/factorial(i));
+end
+printf('value of function at %i is :%f',x0,y_x);
+\ No newline at end of file
diff --git a/1673/CH3/EX3.7/3_7.sce b/1673/CH3/EX3.7/3_7.sce
new file mode 100755
index 000000000..23bef9129
--- /dev/null
+++ b/1673/CH3/EX3.7/3_7.sce
@@ -0,0 +1,9 @@
+//example 3.7
+//interpolation
+//page 89
+clc;clear;close;
+x=[0 1 2 4];
+y=[1 3 9 81];
+//equation is y(5)-4*y(4)+6*y(2)-4*y(2)+y(1)
+y3=(y(4)+6*y(3)-4*y(2)+y(1))/4;
+printf(' the value of missing term of table is :%d',y3);
diff --git a/1673/CH3/EX3.8/3_8.sce b/1673/CH3/EX3.8/3_8.sce
new file mode 100755
index 000000000..0cc288dfa
--- /dev/null
+++ b/1673/CH3/EX3.8/3_8.sce
@@ -0,0 +1,73 @@
+//example 3.8
+//interpolation
+//page 89
+clc;clear;close;
+x=[0.10 0.15 0.20 0.25 0.30];
+y=[0.1003 0.1511 0.2027 0.2553 0.3093];
+h=0.05//interval between values of x
+c=1;
+for i=1:4
+    d1(c)=y(i+1)-y(i);
+    c=c+1;
+end
+c=1;
+for i=1:3
+    d2(c)=d1(i+1)-d1(i);
+    c=c+1
+end
+c=1;
+for i=1:2
+    d3(c)=d2(i+1)-d2(i);
+    c=c+1;
+end
+c=1;
+for i=1:1
+    d4(c)=d3(i+1)-d3(i);
+    c=c+1;
+end
+
+d=[d1(1) d2(1) d3(1) d4(1)];
+x0=0.12;//value at 0.12;
+pp=1;
+y_x=y(1);
+p=(x0-x(1))/h;
+for i=1:4
+    pp=1;
+    for j=1:i
+    pp=pp*(p-(j-1))    
+    end
+y_x=y_x+(pp*d(i))/factorial(i);
+end
+printf('value of function at %f is :%0.4g\n \n',x0,y_x);
+d=[d1(4) d2(3) d3(2) d4(1)];
+x0=0.26;//value at 0.26;
+pp=1;
+y_x=y(5);
+p=(x0-x(5))/h;
+for i=1:4
+    pp=1;
+    for j=1:i
+    pp=pp*(p-(j-1))    
+    end
+y_x=y_x+(pp*d(i))/factorial(i);
+end
+printf('value of function at %f is :%0.4g\n \n',x0,y_x);
+d=[d1(4) d2(3) d3(2) d4(1)];
+x0=0.40;//value at 0.40;
+pp=1;
+y_x=y(5);
+p=(x0-x(5))/h;
+for i=1:4
+    pp=1;
+    for j=1:i
+    pp=pp*(p+(j-1))    
+    end
+y_x=y_x+(pp*d(i))/factorial(i);
+end
+printf('value of function at %f is :%0.4g\n \n',x0,y_x);
+d=[d1(4) d2(3) d3(2) d4(1)];
+x0=0.50;//value at 0.50;
+pp=1;
+y_x=y(5);
+p=(x0-x(5))/h;
+printf('value of function at %f is :%0.5g\n \n',x0,y_x);
diff --git a/1673/CH3/EX3.9/3_9.sce b/1673/CH3/EX3.9/3_9.sce
new file mode 100755
index 000000000..5431b7b70
--- /dev/null
+++ b/1673/CH3/EX3.9/3_9.sce
@@ -0,0 +1,50 @@
+//example 3.9
+//Gauss' forward formula
+//page 3.9
+clc;clear;close;
+x=[1.0 1.05 1.10 1.15 1.20 1.25 1.30];
+y=[2.7183 2.8577 3.0042 3.1582 3.3201 3.4903 3.66693];
+h=0.05//interval between values of x
+c=1;
+for i=1:6
+    d1(c)=y(i+1)-y(i);
+    c=c+1;
+end
+c=1;
+for i=1:5
+    d2(c)=d1(i+1)-d1(i);
+    c=c+1
+end
+c=1;
+for i=1:4
+    d3(c)=d2(i+1)-d2(i);
+    c=c+1;
+end
+c=1;
+for i=1:3
+    d4(c)=d3(i+1)-d3(i);
+    c=c+1;
+end
+c=1;
+for i=1:2
+    d5(c)=d4(i+1)-d4(i);
+    c=c+1;
+end
+c=1;
+for i=1:1
+    d6(c)=d5(i+1)-d5(i);
+    c=c+1;
+end
+d=[d1(4) d2(3) d3(3) d4(2) d5(1) d6(1)];
+x0=1.17;//value at 1.17;
+pp=1;
+y_x=y(4);
+p=(x0-x(4))/h;
+for i=1:6
+    pp=1;
+    for j=1:i
+    pp=pp*(p-(j-1))    
+    end
+y_x=y_x+(pp*d(i))/factorial(i);
+end
+printf('value of function at %f is :%0.4g\n \n',x0,y_x);