From b1f5c3f8d6671b4331cef1dcebdf63b7a43a3a2b Mon Sep 17 00:00:00 2001
From: priyanka
Date: Wed, 24 Jun 2015 15:03:17 +0530
Subject: initial commit / add all books

---
 845/CH7/EX7.1/Ex7_1.sce   | 29 +++++++++++++++++++++++++
 845/CH7/EX7.10/Ex7_10.sce | 23 ++++++++++++++++++++
 845/CH7/EX7.11/Ex7_11.sce | 31 +++++++++++++++++++++++++++
 845/CH7/EX7.12/Ex7_12.sce | 39 ++++++++++++++++++++++++++++++++++
 845/CH7/EX7.13/Ex7_13.sce | 39 ++++++++++++++++++++++++++++++++++
 845/CH7/EX7.14/Ex7_14.sce | 25 ++++++++++++++++++++++
 845/CH7/EX7.15/Ex7_15.sce | 27 ++++++++++++++++++++++++
 845/CH7/EX7.16/Ex7_16.sce | 31 +++++++++++++++++++++++++++
 845/CH7/EX7.2/Ex7_2.sce   | 28 ++++++++++++++++++++++++
 845/CH7/EX7.3/Ex7_3.sce   | 52 +++++++++++++++++++++++++++++++++++++++++++++
 845/CH7/EX7.4/Ex7_4.sce   | 43 +++++++++++++++++++++++++++++++++++++
 845/CH7/EX7.5/Ex7_5.sce   | 28 ++++++++++++++++++++++++
 845/CH7/EX7.6/Ex7_6.sce   | 52 +++++++++++++++++++++++++++++++++++++++++++++
 845/CH7/EX7.7/Ex7_7.sce   | 34 +++++++++++++++++++++++++++++
 845/CH7/EX7.8/Ex7_8.sce   | 54 +++++++++++++++++++++++++++++++++++++++++++++++
 845/CH7/EX7.9/Ex7_9.sce   | 33 +++++++++++++++++++++++++++++
 16 files changed, 568 insertions(+)
 create mode 100755 845/CH7/EX7.1/Ex7_1.sce
 create mode 100755 845/CH7/EX7.10/Ex7_10.sce
 create mode 100755 845/CH7/EX7.11/Ex7_11.sce
 create mode 100755 845/CH7/EX7.12/Ex7_12.sce
 create mode 100755 845/CH7/EX7.13/Ex7_13.sce
 create mode 100755 845/CH7/EX7.14/Ex7_14.sce
 create mode 100755 845/CH7/EX7.15/Ex7_15.sce
 create mode 100755 845/CH7/EX7.16/Ex7_16.sce
 create mode 100755 845/CH7/EX7.2/Ex7_2.sce
 create mode 100755 845/CH7/EX7.3/Ex7_3.sce
 create mode 100755 845/CH7/EX7.4/Ex7_4.sce
 create mode 100755 845/CH7/EX7.5/Ex7_5.sce
 create mode 100755 845/CH7/EX7.6/Ex7_6.sce
 create mode 100755 845/CH7/EX7.7/Ex7_7.sce
 create mode 100755 845/CH7/EX7.8/Ex7_8.sce
 create mode 100755 845/CH7/EX7.9/Ex7_9.sce

(limited to '845/CH7')

diff --git a/845/CH7/EX7.1/Ex7_1.sce b/845/CH7/EX7.1/Ex7_1.sce
new file mode 100755
index 000000000..999247d17
--- /dev/null
+++ b/845/CH7/EX7.1/Ex7_1.sce
@@ -0,0 +1,29 @@
+//Example 7.1
+
+clc
+clear
+
+x = 0:0.2:1;
+y = [1 1.16 3.56 13.96 41.96 101];
+
+n = length(x);
+del = %nan*ones(n,6);
+del(:,1) = y';
+for j = 2:6
+    for i = 1:n-j+1
+        del(i,j) = del(i+1,j-1) - del(i,j-1);
+    end
+end
+del = round(del*10^2)/10^2;
+mprintf("%5s %6s %9s %8s %8s %8s %7s",'x','y','dy','d2y','d3y','d4y','d5y')
+disp([x' del])
+
+h = x(2) - x(1);
+del0 = del(1,:);
+del1 = del(2,:);
+
+df1 = (del1(2) - del1(3)/2 + del1(4)/3 - del1(5)/4) / h;
+d2f0 = (del0(2) - del0(3) + del0(4)*11/12 - del0(5)*5/6) / h^2;
+disp(round(d2f0*10^1)/10^1,"f''''(0) = ")
+disp(round(df1*10)/10,"f''(0.2) = ")
+
diff --git a/845/CH7/EX7.10/Ex7_10.sce b/845/CH7/EX7.10/Ex7_10.sce
new file mode 100755
index 000000000..6c7e359e4
--- /dev/null
+++ b/845/CH7/EX7.10/Ex7_10.sce
@@ -0,0 +1,23 @@
+//Etample 7.10
+
+clc
+clear
+
+t = 0:10:80;
+a = [30 31.63 33.34 35.47 37.75 40.33 43.25 46.69 50.67];
+
+h = t(2) - t(1);
+n = length(t);
+
+Is13 = a(1);
+for i = 2:n-1
+    rem2 = i-fix(i./2).*2;
+    if rem2 == 0 then
+        Is13 = Is13 + 4*a(i);
+    else
+        Is13 = Is13 + 2*a(i);
+    end
+end
+Is13 = (Is13 + a(n))/10^3;
+Is13 = round(h/3*Is13*10^4)/10^4;
+disp(strcat(["v = ",string(Is13)," km/s"]))
diff --git a/845/CH7/EX7.11/Ex7_11.sce b/845/CH7/EX7.11/Ex7_11.sce
new file mode 100755
index 000000000..1ff7ec466
--- /dev/null
+++ b/845/CH7/EX7.11/Ex7_11.sce
@@ -0,0 +1,31 @@
+//Example 7.11
+
+clc
+clear
+
+x = 1:0.1:1.8;
+x = round(x*10)/10;
+y = [1.543 1.669 1.811 1.971 2.151 2.352 2.577 2.828 3.107];
+n = length(x);
+x0 = x(1);
+xn = x(n);
+
+N = [1 2 4 8]
+for j = 1:length(N)
+    h = (xn - x0)./N(j);
+    I = y(1);
+    for xx = x0+h:h:xn-h
+        xx = round(xx*10)/10;
+        I = I + 2*y(x==xx);
+    end
+    Itrap(j) = h/2*(I + y(n));
+    IRomb(1) = Itrap(1);
+    if j~=1 then
+        IRomb(j) = (4^(j-1)*Itrap(j)-IRomb(j-1))/(4^(j-1)-1);
+    end
+end
+IRomb = round(IRomb*10^5)/10^5;
+
+disp(Itrap(length(N)),"Integral using Trapezoidal rule:")
+disp(IRomb(length(N)),"Integral using Romberg''s formula:")
+//In third step of computation of integral using Romberg's formula, author mistakenly took the 1.7672 instead of 1.7684 which resulted in a difference
diff --git a/845/CH7/EX7.12/Ex7_12.sce b/845/CH7/EX7.12/Ex7_12.sce
new file mode 100755
index 000000000..bd73dac92
--- /dev/null
+++ b/845/CH7/EX7.12/Ex7_12.sce
@@ -0,0 +1,39 @@
+//Example 7.12
+
+clc
+clear
+
+function [f] = fun1(x,y)
+    f = 1 / (x+y);
+endfunction
+
+x = 1:0.25:2;
+y = x;
+
+m = length(x);
+n = length(y);
+
+del = %nan*ones(m,n);
+for j = 1:n
+    for i = 1:m
+        del(i,j) = fun1(x(i),y(j));
+    end
+end
+
+hx = x(2) - x(1);
+for i = 1:m
+    I = del(i,1);
+    for j = 2:n-1
+        I = I + 2*del(i,j);
+    end
+    Itrap1(i) = hx/2 * (I+del(i,n));
+end
+Itrap1 = round(Itrap1*10^4)/10^4;
+
+hy = y(2) - y(1);
+Itrap2 = Itrap1(1)
+for i = 2:n-1
+    Itrap2 = Itrap2 + 2* Itrap1(i);
+end
+Itrap2 = round(hy/2*(Itrap2+Itrap1(m))*10^4)/10^4;
+disp(Itrap2,"I = ")
diff --git a/845/CH7/EX7.13/Ex7_13.sce b/845/CH7/EX7.13/Ex7_13.sce
new file mode 100755
index 000000000..bfc8ba495
--- /dev/null
+++ b/845/CH7/EX7.13/Ex7_13.sce
@@ -0,0 +1,39 @@
+//Example 7.13
+
+clc
+clear
+
+function [f] = fun1(x,y)
+    f = sqrt(sin(x+y));
+endfunction
+
+x = 0:%pi/8:%pi/2;
+y = x;
+
+m = length(x);
+n = length(y);
+
+del = %nan*ones(m,n);
+for j = 1:n
+    for i = 1:m
+        del(i,j) = fun1(x(i),y(j));
+    end
+end
+
+hx = x(2) - x(1);
+for i = 1:m
+    I = del(i,1);
+    for j = 2:n-1
+        I = I + 2*del(i,j);
+    end
+    Itrap1(i) = hx/2 * (I+del(i,n));
+end
+Itrap1 = round(Itrap1*10^4)/10^4;
+
+hy = y(2) - y(1);
+Itrap2 = Itrap1(1)
+for i = 2:n-1
+    Itrap2 = Itrap2 + 2* Itrap1(i);
+end
+Itrap2 = round(hy/2*(Itrap2+Itrap1(m))*10^4)/10^4;
+disp(Itrap2,"I = ")
diff --git a/845/CH7/EX7.14/Ex7_14.sce b/845/CH7/EX7.14/Ex7_14.sce
new file mode 100755
index 000000000..47da1f96f
--- /dev/null
+++ b/845/CH7/EX7.14/Ex7_14.sce
@@ -0,0 +1,25 @@
+//Example 7.14
+
+clc
+clear
+
+n = 1;
+if n==1 then
+    M = [0 2];
+elseif n==2
+    M = [sqrt(1/3) 1; -sqrt(1/3) 1];
+elseif n==3
+    M = [0 8/9; -0.774597 5/9; 0.774597 5/9];
+elseif n==4
+    M = [-0.339981 0.652145; -0.861136 0.347855; 0.339981 0.652145; 0.861136 0.347855];
+elseif n==5
+    M = [-0 0.568889; -0.538469 0.467914; -0.906180 0.236927; 0 0.568889; 0.538469 0.467914; 0.906180 0.236927];
+elseif n==6
+    M = [-0.238619 0.467914; -0.661209 0.360762; -0.932470 0.171325; 0.238619 0.467914; 0.661209 0.360762; 0.932470 0.171325];
+end
+
+X = M(:,1);
+W = M(:,2);
+
+disp(X,"E1 = ")
+disp(W,"W1 = ")
diff --git a/845/CH7/EX7.15/Ex7_15.sce b/845/CH7/EX7.15/Ex7_15.sce
new file mode 100755
index 000000000..b3ec6213b
--- /dev/null
+++ b/845/CH7/EX7.15/Ex7_15.sce
@@ -0,0 +1,27 @@
+//Example 7.15
+
+clc
+clear
+
+n = 2;
+if n==1 then
+    M = [0 2];
+elseif n==2
+    M = [sqrt(1/3) 1; -sqrt(1/3) 1];
+elseif n==3
+    M = [0 8/9; -0.774597 5/9; 0.774597 5/9];
+elseif n==4
+    M = [-0.339981 0.652145; -0.861136 0.347855; 0.339981 0.652145; 0.861136 0.347855];
+elseif n==5
+    M = [-0 0.568889; -0.538469 0.467914; -0.906180 0.236927; 0 0.568889; 0.538469 0.467914; 0.906180 0.236927];
+elseif n==6
+    M = [-0.238619 0.467914; -0.661209 0.360762; -0.932470 0.171325; 0.238619 0.467914; 0.661209 0.360762; 0.932470 0.171325];
+end
+
+X = M(:,1);
+W = M(:,2);
+
+disp(W(1),"W1 = ")
+disp(W(2),"W2 = ")
+disp(X(1),"E1 = ")
+disp(X(2),"E2 = ")
diff --git a/845/CH7/EX7.16/Ex7_16.sce b/845/CH7/EX7.16/Ex7_16.sce
new file mode 100755
index 000000000..fcd56815d
--- /dev/null
+++ b/845/CH7/EX7.16/Ex7_16.sce
@@ -0,0 +1,31 @@
+//Example 7.16
+
+clc
+clear
+
+function [f] = fun1(x)
+    f = 3*x^2 + x^3;
+endfunction
+
+n = 4;
+if n==1 then
+    M = [0 2];
+elseif n==2
+    M = [sqrt(1/3) 1];
+elseif n==3
+    M = [0 8/9; -0.774597 5/9; 0.774597 5/9];
+elseif n==4
+    M = [-0.339981 0.652145; -0.861136 0.347855; 0.339981 0.652145; 0.861136 0.347855];
+elseif n==5
+    M = [-0 0.568889; -0.538469 0.467914; -0.906180 0.236927; 0 0.568889; 0.538469 0.467914; 0.906180 0.236927];
+elseif n==6
+    M = [-0.238619 0.467914; -0.661209 0.360762; -0.932470 0.171325; 0.238619 0.467914; 0.661209 0.360762; 0.932470 0.171325];
+end
+
+X = M(:,1);
+W = M(:,2);
+I = 0;
+for i = 1:length(X)
+    I = I + W(i)*fun1(X(i));
+end
+disp(I,"I = ")
diff --git a/845/CH7/EX7.2/Ex7_2.sce b/845/CH7/EX7.2/Ex7_2.sce
new file mode 100755
index 000000000..01a902a7d
--- /dev/null
+++ b/845/CH7/EX7.2/Ex7_2.sce
@@ -0,0 +1,28 @@
+//Example 7.2
+
+clc
+clear
+
+x = 1.4:0.2:2.2;
+y = [4.0552 4.953 6.0496 7.3891 9.025];
+
+n = length(x);
+del = %nan*ones(n,5);
+del(:,1) = y';
+for j = 2:5
+    for i = 1:n-j+1
+        del(i+j-1,j) = del(i+j-1,j-1) - del(i+j-2,j-1);
+    end
+end
+mprintf("%5s %6s %10s %10s %9s %9s",'x','y','dy','d2y','d3y','d4y')
+disp([x' del])
+
+h = x(2) - x(1);
+deln = del(5,:);
+
+dfn = (deln(2) + deln(3)/2 + deln(4)/3 + deln(5)/4) / h;
+d2fn = (deln(3) + deln(4) + deln(5)*11/12) / h^2;
+dfn = round(dfn*10^4)/10^4;
+d2fn = round(d2fn*10^4)/10^4;
+disp(dfn,"y''(2.2) = ")
+disp(d2fn,"y''''(2.2) = ")
diff --git a/845/CH7/EX7.3/Ex7_3.sce b/845/CH7/EX7.3/Ex7_3.sce
new file mode 100755
index 000000000..1d02509df
--- /dev/null
+++ b/845/CH7/EX7.3/Ex7_3.sce
@@ -0,0 +1,52 @@
+//Example 7.3
+
+clc
+clear
+
+x = 0:4;
+y = [6.9897 7.4036 7.7815 8.1281 8.451];
+
+n = length(x);
+del = %nan*ones(n,5);
+del(:,1) = y';
+for j = 2:6
+    for i = 1:n-j+1
+        del(i,j) = del(i+1,j-1) - del(i,j-1);
+    end
+end
+del(:,1) = [];
+n0 = length(del(1,:));
+
+X = 2;
+i = find(x==X);
+dowy = 0;
+
+for j = 1:n0
+    if j==2*int(j/2) then
+        add = del(i,j);
+    else
+        add = (del(i-1,j) + del(i,j))/2;
+        i = i-1;
+        if i==0 then
+            break
+        end
+    end
+
+    if add == %nan then
+        break
+    else
+        dowy(j) = add;
+    end
+end
+mprintf("%5s %6s %10s %9s %9s %9s",'x','y','dy','d2y','d3y','d4y')
+disp([x' y' del])
+
+mu = 1;
+h = x(2) - x(1);
+dy2 = mu/h*(dowy(1) - 1/6*dowy(3));
+d2y2 = mu/h^2*(dowy(2)-1/12*dowy(4));
+dy2 = round(dy2*10^4)/10^4;
+d2y2 = round(d2y2*10^4)/10^4;
+
+disp(dy2,"y''(2) = ")
+disp(d2y2,"y''''(2) = ")
diff --git a/845/CH7/EX7.4/Ex7_4.sce b/845/CH7/EX7.4/Ex7_4.sce
new file mode 100755
index 000000000..f35e51791
--- /dev/null
+++ b/845/CH7/EX7.4/Ex7_4.sce
@@ -0,0 +1,43 @@
+//Example 7.4
+
+clc
+clear
+
+x = [0.15 0.21 0.23 0.27 0.32 0.35];
+y = [0.1761 0.3222 0.3617 0.4314 0.5051 0.5441];
+
+n = length(x);
+del = %nan*ones(n,6);
+del(:,1) = y';
+for j = 2:6
+    for i = 1:n-j+1
+        del(i,j) = (del(i+1,j-1) - del(i,j-1)) / (x(i+j-1)-x(i));
+    end
+end
+del(:,1) = [];
+del = round(del*10^3)/10^3;
+mprintf("%5s %6s %10s %10s %8s %9s %9s",'x','y','dy','d2y','d3y','d4y','d5y')
+disp([x' y' del])
+X = poly(0, "X");
+del0 = del(1,:);
+y0 = y(1);
+Y = y0;
+for i = 1:length(del0)
+    p = 1;
+    for j = 1:i
+        p = p*(X-x(j));
+    end
+    Y = Y + p*del0(i);
+end
+
+dydx = derivat(Y);
+d2ydx = derivat(dydx);
+
+XX = 0.25;
+dy = horner(dydx,XX);
+d2y = horner(d2ydx,XX);
+
+disp(round(dy*10^4)/10^4,"y''(0.25) = ")
+disp(d2ydx,"y''''(x) = ")
+disp(d2y,"y''''(0.25) = ")
+//The constant term in y''(x) is incorrectly computed to -91.7 instead of -97.42 in the text.
diff --git a/845/CH7/EX7.5/Ex7_5.sce b/845/CH7/EX7.5/Ex7_5.sce
new file mode 100755
index 000000000..ebc95f32b
--- /dev/null
+++ b/845/CH7/EX7.5/Ex7_5.sce
@@ -0,0 +1,28 @@
+//Example 7.5
+
+clc
+clear
+
+function [f] = y(x)
+    f = -1/x;
+endfunction
+
+H = [0.0128 0.0064 0.0032];
+n = length(H);
+x = 0.05;
+h = H(1);
+Fh = (y(x+h) - y(x-h)) / (2*h);
+Fh2 = (y(x+h/2) - y(x-h/2)) / (h);
+Fh4 = (y(x+h/4) - y(x-h/4)) / (h/2);
+
+F1h2 = (4*Fh2 - Fh) / (4-1);
+F1h4 = (4*Fh4 - Fh2) / (4-1);
+F2h4 = (4^2*F1h4 - F1h2) / (4^2-1);
+del = %nan*ones(n,3);
+del(:,1) = [Fh Fh2 Fh4]';
+del(1:2,2) = [F1h2 F1h4]';
+del(1,3) = F2h4;
+
+disp(del(1,n),"y''(0.05) = ")
+Exact = 1/x^2;
+disp(Exact,"Exact Value:")
diff --git a/845/CH7/EX7.6/Ex7_6.sce b/845/CH7/EX7.6/Ex7_6.sce
new file mode 100755
index 000000000..d8b00fe69
--- /dev/null
+++ b/845/CH7/EX7.6/Ex7_6.sce
@@ -0,0 +1,52 @@
+//Example 7.6
+
+clc
+clear
+
+function [I] = trap (fun,a,b,n)
+// Integrate the function over the interval using Trapezoidal Formula
+// trap (fun,a,b,n)
+// fun - function to be integrated
+// a - lower limit of integration
+// b - upper limit of integration
+// n - No. of times trapezoidal rule needs to be performed
+
+N = n + 1;  // N - total no. of points
+h = (b-a) / (N-1);
+x = linspace(a,b,N);
+y = fun(x);
+
+sum1 = y(1) + 2 * sum(y(2:N-1)) + y(N);
+I = h * sum1 / 2;                   // Trapezoidal Integral Value
+endfunction
+
+function [I] = simp13 (fun,a,b,n)
+// Integrate the function over the interval using Simpson's 1/3rd rule
+// simp13 (fun,a,b,n)
+// fun - function to be integrated
+// a - lower limit of integration
+// b - upper limit of integration
+// n - No. of times simpson's 1/3rd rule needs to be performed
+
+N = 2 * n + 1;      // N - total no. of points
+h = (b-a) / (N-1);
+x = linspace(a,b,N);
+y = fun(x);
+
+sum1 = y(1) + 4 * sum(y(2:2:N-1)) + 2 * sum(y(3:2:N-2)) + y(N);
+I = h* sum1 / 3;                    // Simpson's 1/3rd Integral Value
+endfunction
+
+n = 6;
+ntrap = n;
+ns13 = n/2;
+I = [trap(sin,0,%pi,ntrap); simp13(sin,0,%pi,ns13)];
+I = round(I*10^4)/10^4;
+true = integrate('sin(x)','x',0,%pi);
+err = abs(true - I) / true*100;
+err = round(err*100)/100;
+
+disp(I(1),"y_trap = ")
+disp(I(2),"y_simp13 = ")
+disp(err(1),"error_trap = ")
+disp(err(2),"error_simp13 = ")
diff --git a/845/CH7/EX7.7/Ex7_7.sce b/845/CH7/EX7.7/Ex7_7.sce
new file mode 100755
index 000000000..df21dadbd
--- /dev/null
+++ b/845/CH7/EX7.7/Ex7_7.sce
@@ -0,0 +1,34 @@
+//Example 7.7
+
+clc
+clear
+
+function [I] = simp13 (fun,a,b,n)
+// Integrate the function over the interval using Simpson's 1/3rd rule
+// simp13 (fun,a,b,n)
+// fun - function to be integrated
+// a - lower limit of integration
+// b - upper limit of integration
+// n - No. of times simpson's 1/3rd rule needs to be performed
+
+N = 2 * n + 1;      // N - total no. of points
+h = (b-a) / (N-1);
+x = linspace(a,b,N);
+y = fun(x);
+
+sum1 = y(1) + 4 * sum(y(2:2:N-1)) + 2 * sum(y(3:2:N-2)) + y(N);
+I = h* sum1 / 3;                    // Simpson's 1/3rd Integral Value
+endfunction
+
+n = 8;
+ns13 = n/2;
+I = simp13(log,1,5,ns13);
+I = round(I*10^4)/10^4;
+deff('[y] = true(x)',['y = x * log(x) - x']);
+trueVal = true(5) - true(1);
+err = abs(trueVal - I) / trueVal*100;
+err = round(err*100)/100;
+
+disp(I,"y_simp13 = ")
+disp(trueVal,"Actual Integral = ")
+disp(err,"error_simp13 = ")
diff --git a/845/CH7/EX7.8/Ex7_8.sce b/845/CH7/EX7.8/Ex7_8.sce
new file mode 100755
index 000000000..0318b5ee7
--- /dev/null
+++ b/845/CH7/EX7.8/Ex7_8.sce
@@ -0,0 +1,54 @@
+//Example 7.8
+
+clc
+clear
+
+function [I] = trap (fun,a,b,n)
+// Integrate the function over the interval using Trapezoidal Formula
+// trap (fun,a,b,n)
+// fun - function to be integrated
+// a - lower limit of integration
+// b - upper limit of integration
+// n - No. of times trapezoidal rule needs to be performed
+
+N = n + 1;  // N - total no. of points
+h = (b-a) / (N-1);
+x = linspace(a,b,N);
+y = fun(x);
+
+sum1 = y(1) + 2 * sum(y(2:N-1)) + y(N);
+I = h * sum1 / 2;                   // Trapezoidal Integral Value
+endfunction
+
+function [I] = simp13 (fun,a,b,n)
+// Integrate the function over the interval using Simpson's 1/3rd rule
+// simp13 (fun,a,b,n)
+// fun - function to be integrated
+// a - lower limit of integration
+// b - upper limit of integration
+// n - No. of times simpson's 1/3rd rule needs to be performed
+
+N = 2 * n + 1;      // N - total no. of points
+h = (b-a) / (N-1);
+x = linspace(a,b,N);
+y = fun(x);
+
+sum1 = y(1) + 4 * sum(y(2:2:N-1)) + 2 * sum(y(3:2:N-2)) + y(N);
+I = h* sum1 / 3;                    // Simpson's 1/3rd Integral Value
+endfunction
+
+function [f] = fun1(x)
+    f = 1 ./(1+x^2);
+endfunction
+
+
+n = 4;
+ntrap = n;
+ns13 = n/2;
+I = [trap(fun1,0,1,ntrap); simp13(fun1,0,1,ns13)];
+I = round(I*10^4)/10^4;
+true = intg(0,1,fun1);
+
+disp(I(1),"y_trap = ")
+disp(I(2),"y_simp13 = ")
+disp(I(2)*4,"Approx pi = ")
diff --git a/845/CH7/EX7.9/Ex7_9.sce b/845/CH7/EX7.9/Ex7_9.sce
new file mode 100755
index 000000000..ad4ef22ae
--- /dev/null
+++ b/845/CH7/EX7.9/Ex7_9.sce
@@ -0,0 +1,33 @@
+//Example 7.9
+
+clc
+clear
+
+function [I] = simp13 (fun,a,b,n)
+// Integrate the function over the interval using Simpson's 1/3rd rule
+// simp13 (fun,a,b,n)
+// fun - function to be integrated
+// a - lower limit of integration
+// b - upper limit of integration
+// n - No. of times simpson's 1/3rd rule needs to be performed
+
+N = 2 * n + 1;      // N - total no. of points
+h = (b-a) / (N-1);
+x = linspace(a,b,N);
+y = fun(x);
+
+sum1 = y(1) + 4 * sum(y(2:2:N-1)) + 2 * sum(y(3:2:N-2)) + y(N);
+I = h* sum1 / 3;                    // Simpson's 1/3rd Integral Value
+endfunction
+
+function [f] = fun1(x)
+    f = sqrt(2/%pi)*exp(-x^2/2);
+endfunction
+
+h = 0.125;
+n = (1-0)/h;
+ns13 = n/2;
+I = simp13(fun1,0,1,ns13);
+I = round(I*10^4)/10^4;
+
+disp(I,"Integral value, I = ")
-- 
cgit