19 files changed, 883 insertions, 0 deletions
diff --git a/260/CH13/EX13.1/13_1.sce b/260/CH13/EX13.1/13_1.sce
new file mode 100644
index 000000000..848a82c15
--- /dev/null
+++ b/260/CH13/EX13.1/13_1.sce
@@ -0,0 +1,20 @@
+//Eg-13.1
+//pg-522
+
+clear
+clc
+
+//Given equation dy/dx = x. Hence F'(x) = 1 and F''(x) = 0 
+
+//In the taylor series expnsion, if we write y0 in place of y(x0), we have
+// y(1) = y1 = y(x0+h) = y0 + hx0 + h^2/2
+
+x0 = 0;
+y0 = 0;    //Initial condition
+
+h = 1;    //taking the value ourself
+
+y1 = y0 + h*x0 + h^2/2;
+
+printf('The value of y at x = 1 is %f\n',y1)
+printf(' This is the exact solution since the higher derivatives starting from second order derivative of F vanish\n')
+\ No newline at end of file
diff --git a/260/CH13/EX13.10/13_10.sce b/260/CH13/EX13.10/13_10.sce
new file mode 100644
index 000000000..f789f08fb
--- /dev/null
+++ b/260/CH13/EX13.10/13_10.sce
@@ -0,0 +1,32 @@
+//find the function like eig in matlab
+//Eg-13.10
+//pg-548
+
+clear
+clc
+
+//From the given differential equations we have,
+//    F0 = y0 - y1 + 2*y2^2;
+//    F1 = 3*y0 + y1^2 - y2;
+//    F2 = 2*y0 + 0.5*y1^2 - 0.5*y2^2;
+
+//Therefore, the jacobian matrix is 
+
+J = [1 -1 4;3 2 -1;2 1 -1]
+Ji = inv(J)
+
+printf('J = ')
+
+disp(J)
+
+printf('The largest eigenvalue can be computed by the power method(see chapter 5). It is found \nto be 3. To find the minimum eigen value, we determine the inverse of this matrix\n')
+printf('\nInverse of J, Ji = \n')
+
+disp(Ji)
+
+//ro = Lmax/Lmin 
+
+ro = 3/1;
+
+printf('\nThe power method gives the largest eigenvalue of this matrix to be 1. Therefore, the \nminumum eigenvalue of J is 1/1 = 1. The stiffness ratio, therefore is,\n\n    ro = %f\n',ro)
+
diff --git a/260/CH13/EX13.11/13_11.sce b/260/CH13/EX13.11/13_11.sce
new file mode 100644
index 000000000..ee3732114
--- /dev/null
+++ b/260/CH13/EX13.11/13_11.sce
@@ -0,0 +1,61 @@
+//Eg-13.11
+//pg-550
+
+clear
+clc
+
+deff('out = func(in1,in2)','out = -2*in2')
+
+for(j = 1:3)
+    
+    x(1) = 0;
+    y(1,j) = 1;
+    yb(1,j) = 0;
+    
+end
+
+
+h(1) = 0.43;
+h(2) = 0.3;
+h(3) = 0.03;
+
+z(1) = exp(-2*x(1));
+
+for(i = 2:11)
+    
+    x(i) = 0.1+(i-2)*0.1;
+    z(i+1) = exp(-2*x(i));
+end
+
+for(j = 1:3)
+    
+    
+    for(i = 1:10)
+        yb(i+1,j) = y(i) + h(j)*func(x(i),y(i,j));
+        y(i+1,j) = y(i) + h(j)/2*(func(x(i),y(i,j)) + func(x(i+1),yb(i+1,j)));
+    end
+
+
+end
+
+printf('   x           h1         h2          h3         exact\n')
+
+for(i = 2:11)
+    printf('%f    %f    %f    %f    %f\n',x(i),y(i,1),y(i,2),y(i,3),z(i+1))
+end
+
+for(i = 1:10)
+    a(i) = x(i+1);
+    b(i) = y(i+1,1);
+    c(i) = y(i+1,2);
+    d(i) = y(i+1,3);
+    e(i) = z(i+2);
+end
+
+clf()
+//plot(a,[b c d e])
+plot(a,b,'s')
+plot(a,c,'o')
+plot(a,d,'d')
+plot(a,e,'-')
+legend('h = 0.43','h = 0.3','h = 0.03','analytical')
+\ No newline at end of file
diff --git a/260/CH13/EX13.12/13_12.sce b/260/CH13/EX13.12/13_12.sce
new file mode 100644
index 000000000..4659f5353
--- /dev/null
+++ b/260/CH13/EX13.12/13_12.sce
@@ -0,0 +1,32 @@
+//Eg-13.12
+//pg-552
+
+clear
+clc
+
+y(1) = 1;
+
+h = 1;
+
+printf('For h = 1\n')
+printf('              x           y\n')
+
+for(i = 1:3)
+    x(i) = i;    //since h = 1
+    y(i+1) = y(i)/(1+2*h);
+    printf('            %f    %f\n',x(i),y(i+1))
+end
+
+h = 0.5;
+printf('\nFor h = 0.5\n')
+printf('              x           y\n')
+
+n = (3.0-0.5)/0.5+1;
+
+for(i = 1:n)
+    x(i) = 0.5 + (i-1)*h;    //since h = 0.5
+    y(i+1) = y(i)/(1+2*h);
+    printf('            %f    %f\n',x(i),y(i+1))
+end
+
+printf('Observe that the implicit method is stable for h = 1, whereas the explicit method is not.')
+\ No newline at end of file
diff --git a/260/CH13/EX13.13/13_13.sce b/260/CH13/EX13.13/13_13.sce
new file mode 100644
index 000000000..6d2c5dff4
--- /dev/null
+++ b/260/CH13/EX13.13/13_13.sce
@@ -0,0 +1,50 @@
+//Eg-13.13
+//pg-554
+
+clear
+clc
+
+y(1) = 1;
+
+deff('out = func(in1,in2)','out = 1-in1+4*in2')
+
+h = 0.1;
+
+//The index again is 1-11 instead of 0-10
+
+for(i = 1:11)
+    
+    x(i) = 0 + (i-1)*h;
+    
+end
+
+a = (2^0.5-1)/2;
+b = (2-2^0.5)/2;
+c = -(2^0.5)/2;
+d = 1 + (2^0.5)/2;
+
+
+for(i = 1)
+    k1(i) = h*func(x(i),y(i));
+    k2(i) = h*func(x(i)+h/2,y(i)+k1(i)/2);
+    k3(i) = h*func(x(i)+h/2,y(i)+a*k1(i)+b*k2(i));
+    k4(i) = h*func(x(i)+h,y(i)+c*k2(i)+d*k3(i));
+    y(i+1) = y(i) + 1/6*(k1(i)+2*b*k2(i)+2*d*k3(i)+k4(i));
+end
+
+printf('The solution is summarized as :\n')
+printf('   x            y\n')
+for(i = 1:2)
+    
+    printf('%f    %f\n',x(i),y(i))
+
+end
+//we take x(1) = 0 and x(2) = 0.1, y(1) = 1 and y(2) = 1.608933
+
+h = 0.1;
+
+for(i = 3:5)
+    yb(i) = y(i-2) + 2*0.1*func(x(i-1),y(i-1));
+    y(i) = y(i-1) + 0.1/2*(func(x(i-1),y(i-1)) + func(x(i),yb(i)))
+    printf('%f    %f\n',x(i),y(i))
+end
diff --git a/260/CH13/EX13.14/13_14.sce b/260/CH13/EX13.14/13_14.sce
new file mode 100644
index 000000000..7c857b4a6
--- /dev/null
+++ b/260/CH13/EX13.14/13_14.sce
@@ -0,0 +1,40 @@
+//Eg-13.14
+//pg-556
+
+clear
+clc
+
+deff('out = func(in1,in2)','out = -0.3*(in2-50)^1.25')
+
+h = 0.1;
+x(1) = 0;
+y(1) = 100;
+y(2) = 96.20249;    //from the question
+x(2) = x(1) + h;
+
+n = 10/h;
+
+for(i = 1:n)
+    
+    x(i+1) = x(i) + h ;
+    
+end
+
+
+
+
+for(i = 2:n)
+
+    yb(i+1) = y(i-1) + 2*h*func(x(i),y(i));
+    y(i+1) = y(i) + h/2*(func(x(i),y(i)) + func(x(i+1),yb(i+1)));
+    
+end
+
+printf('  x            y\n')
+for(i = 1:10)
+   
+    printf('%f    %f\n',x(i*10+1),y(i*10+1))
+
+end
+
+printf('\n\nThe value of T at t = 10 agrees with that of the analytical solution:\nT = 50 + (0.075t + 0.37604)^-4\n')
+\ No newline at end of file
diff --git a/260/CH13/EX13.15/13_15.sce b/260/CH13/EX13.15/13_15.sce
new file mode 100644
index 000000000..b511cc26b
--- /dev/null
+++ b/260/CH13/EX13.15/13_15.sce
@@ -0,0 +1,59 @@
+//Eg-13.15
+//pg-559
+
+clear
+clc
+
+deff('out = func(in1,in2)','out = 0.5*(1+in1)*in2^2')
+
+x(1) = 0;
+y(1) = 1;
+
+h = 0.001;
+
+n = 0.3/0.001;
+
+for(i = 1:n)
+    x(i+1) = x(i) + h;
+end
+
+a = (2^0.5-1)/2;
+b = (2-2^0.5)/2;
+c = -(2^0.5)/2;
+d = 1 + (2^0.5)/2;
+
+for(i = 1:n)
+    
+    k1(i) = h*func(x(i),y(i));
+    k2(i) = h*func(x(i)+h/2,y(i)+k1(i)/2);
+    k3(i) = h*func(x(i)+h/2,y(i)+a*k1(i)+b*k2(i));
+    k4(i) = h*func(x(i)+h,y(i)+c*k2(i)+d*k3(i));
+    y(i+1) = y(i) + 1/6*(k1(i)+2*b*k2(i)+2*d*k3(i)+k4(i));
+    
+    
+end
+
+printf('   t          y          F(t,y)\n')
+
+for(i = 0:3)
+    printf('%f    %f    %f\n',x(i*100+1),y(i*100+1),func(x(i*100+1),y(i*100+1)))
+end
+
+//Using equations [61] and [62]
+
+y0 = y(1);
+y1 = y(101);
+y2 = y(201);
+y3 = y(301);
+x0 = x(1);
+x1 = x(101);
+x2 = x(201);
+x3 = x(301);
+x4 = 0.4;
+h = 0.1;
+
+y4b = y0 + 4*h/3*(2*func(x3,y3) - func(x2,y2) + 2*func(x1,y1));
+
+y4 = y2 + h/3*(func(x2,y2) + 4*func(x3,y3) + func(x4,y4b));
+
+printf('\nThe value of y4 = %f\n',y4)
+\ No newline at end of file
diff --git a/260/CH13/EX13.16/13_16.sce b/260/CH13/EX13.16/13_16.sce
new file mode 100644
index 000000000..23da481da
--- /dev/null
+++ b/260/CH13/EX13.16/13_16.sce
@@ -0,0 +1,52 @@
+//Eg-13.16
+//pg-562
+
+clear
+clc
+
+deff('out = func(in1,in2)','out = sin(in1)-in2')
+
+//To get the answers in the text book we shouldn't take sint=Vs to be unity!
+
+x(1) = 0;
+y(1) = 0;
+
+h = 0.1;
+
+for(i = 1:10)
+    x(i+1) = x(i) + h;
+end
+
+a = (2^0.5-1)/2;
+b = (2-2^0.5)/2;
+c = -(2^0.5)/2;
+d = 1 + (2^0.5)/2;
+
+printf('   t          y\n')
+
+for(i = 1:4)
+    
+    k1(i) = h*func(x(i),y(i));
+    k2(i) = h*func(x(i)+h/2,y(i)+k1(i)/2);
+    k3(i) = h*func(x(i)+h/2,y(i)+a*k1(i)+b*k2(i));
+    k4(i) = h*func(x(i)+h,y(i)+c*k2(i)+d*k3(i));
+    y(i+1) = y(i) + 1/6*(k1(i)+2*b*k2(i)+2*d*k3(i)+k4(i));
+    printf('%f    %f\n',x(i),y(i))
+    
+end
+
+h = 0.1;
+
+for(i = 4:10)
+    yb(i+1) = y(i-3) + 4*h/3*(2*func(x(i-1),y(i-1)) + 2*func(x(i-2),y(i-2)));
+    y(i+1) = y(i-1) + h/3*(func(x(i-1),y(i-1)) + 4*func(x(i),y(i)) + func(x(i+1),yb(i+1)) );
+
+    
+end
+
+printf('\n\nUsing the Milnes predictor-corrector Method:\n\n')
+
+printf('   t          y\n')
+for(i = 5:11)
+    printf('%f    %f\n',x(i),y(i))
+end
diff --git a/260/CH13/EX13.17/13_17.sce b/260/CH13/EX13.17/13_17.sce
new file mode 100644
index 000000000..45424b296
--- /dev/null
+++ b/260/CH13/EX13.17/13_17.sce
@@ -0,0 +1,56 @@
+//Eg-13.17
+//pg-563
+
+clear
+clc
+
+deff('out = func(in1,in2)','out = in1^3*(cos(in2))^2 - in1*sin(2*in2)')
+// Using the analytical expression in the question
+deff('out = F(in1)','out = atan(0.5*(in1^2 + exp(-in1^2) - 1))')
+
+
+x(1) = 0;
+y(1) = 0;
+
+h = 0.1;
+
+for(i = 1:10)
+    x(i+1) = x(i) + h;
+    
+end
+
+//Using RK method
+
+a = (2^0.5-1)/2;
+b = (2-2^0.5)/2;
+c = -(2^0.5)/2;
+d = 1 + (2^0.5)/2;
+
+printf('   t          y         F(x,y)\n')
+
+for(i = 1:4)
+    
+    k1(i) = h*func(x(i),y(i));
+    k2(i) = h*func(x(i)+h/2,y(i)+k1(i)/2);
+    k3(i) = h*func(x(i)+h/2,y(i)+a*k1(i)+b*k2(i));
+    k4(i) = h*func(x(i)+h,y(i)+c*k2(i)+d*k3(i));
+    y(i+1) = y(i) + 1/6*(k1(i)+2*b*k2(i)+2*d*k3(i)+k4(i));
+    printf('%f    %f    %f\n',x(i),y(i),func(x(i),y(i)))
+    
+end
+
+printf('\n\nUsing Adams Predictor-Corrector Method & its comparision with the analytical solution:\n\n')
+
+h = 0.1;
+for(i = 4:6)
+    
+    yb(i+1) = y(i) + h/24*(55*func(x(i),y(i)) - 59*func(x(i-1),y(i-1)) + 37*func(x(i-2),y(i-2)) - 9*func(x(i-3),y(i-3)));
+    y(i+1) = y(i) + h/24*(9*func(x(i+1),yb(i+1)) + 19*func(x(i),y(i)) - 5*func(x(i-1),y(i-1)) + func(x(i-2),y(i-2)));
+    
+end
+
+printf('   t         yadams    yanalytical\n')
+
+for(i = 5:6)
+    printf('%f    %f    %f\n',x(i),y(i),F(x(i)))
+end
diff --git a/260/CH13/EX13.18/13_18.sce b/260/CH13/EX13.18/13_18.sce
new file mode 100644
index 000000000..096181d62
--- /dev/null
+++ b/260/CH13/EX13.18/13_18.sce
@@ -0,0 +1,61 @@
+//Eg-13.18
+//pg-565
+
+clear
+clc
+
+deff('out = func(in1,in2)','out = (1-in2*sin(in1)/cos(in1))')
+
+x(1) = 0;
+y(1) = 0;
+
+h = 0.1;
+
+for(i = 1:10)
+    x(i+1) = x(i) + h;
+    
+end
+
+//Using RK method
+
+a = (2^0.5-1)/2;
+b = (2-2^0.5)/2;
+c = -(2^0.5)/2;
+d = 1 + (2^0.5)/2;
+
+printf('   x          y\n')
+
+for(i = 1:4)
+    
+    k1(i) = h*func(x(i),y(i));
+    k2(i) = h*func(x(i)+h/2,y(i)+k1(i)/2);
+    k3(i) = h*func(x(i)+h/2,y(i)+a*k1(i)+b*k2(i));
+    k4(i) = h*func(x(i)+h,y(i)+c*k2(i)+d*k3(i));
+    y(i+1) = y(i) + 1/6*(k1(i)+2*b*k2(i)+2*d*k3(i)+k4(i));
+    
+    
+end
+
+for(i = 2:4)
+    printf('%f    %f\n',x(i),y(i))
+end
+
+x4 = x(5);
+x3 = x(4);
+x2 = x(3);
+x1 = x(2);
+x0 = x(1);
+y3 = y(4);
+y2 = y(3);
+y1 = y(2);
+y0 = y(1);
+
+//Using the equation : y4 = 1/25*[48*y3 - 36*y2 + 16*y1 - 3*y0 + 12*h*func(x4,y4)]
+// => y4 - 1/25*12*h*func(x4,y4) = 1/25*[48*y3 - 36*y2 + 16*y1 - 3*y0 ]
+// => y4 = 0.39731/1.02029    analytically
+
+y4 = 0.39731/1.02029;
+
+printf('\n\nThe value of y4 = %f\n',y4)
+
+printf('\nThe analytical solution of this problem is y = sinx.\nThus, the exact solution is y(x = 0.4) = %f\n',y4)
+\ No newline at end of file
diff --git a/260/CH13/EX13.19/13_19.sce b/260/CH13/EX13.19/13_19.sce
new file mode 100644
index 000000000..918639208
--- /dev/null
+++ b/260/CH13/EX13.19/13_19.sce
@@ -0,0 +1,42 @@
+//Eg-13.19
+//pg-567
+
+clear
+clc
+
+deff('out = func(in1,in2,in3)','out = in3 - in2')
+
+//At t = 0
+
+t = 0;
+y0(1) = 1;
+y1(1) = y0 - 1;
+
+//Now all the values are known at t = 0. We will now use the Euler's method to compute y0 at (t+h), with h = 0.1
+
+//At t = 0.1 
+
+t = 0.1;
+y0(2) = y0(1) + t*func(0,y0(1),y1(1));
+y1(2) = y0(2) - 1;
+
+//Similarly
+
+t = 0.2;
+y0(3) = y0(2) + t*func(0.1,y0(2),y1(2));
+y1(3) = y0(3) - 1;
+
+t = 0.3;
+y0(4) = y0(3) + t*func(0.2,y0(3),y1(3));
+y1(4) = y0(4) - 1;
+
+printf('Therefore at t = 0.3, y0 = %f, y1 = %f\n',y0(3),y1(3))
+
+//printf('%f    %f\n%f    %f\n%f    %f\n',y0(1),y1(1),y0(2),y1(2),y0(3),y1(3))
+
+
+
+
+
+
+
diff --git a/260/CH13/EX13.2/13_2.sce b/260/CH13/EX13.2/13_2.sce
new file mode 100644
index 000000000..06af7e065
--- /dev/null
+++ b/260/CH13/EX13.2/13_2.sce
@@ -0,0 +1,26 @@
+//Eg-13.2
+//pg-523
+
+clc
+clear
+
+
+//F = x^2*y = dy/dx
+
+deff('[out] = func(in1,in2)','out = in1^2*in2')
+
+y(1) = 1;    //Initial condition
+x(1) = 0;
+z(1) = exp(x(1)^3/3);
+
+h = 0.1;
+
+printf('x           yEuler      yexact\n')
+for(i = 1:10)
+    x(i) = 0.1*i;
+    y(i+1) = y(i) + h*func(x(i),y(i));
+    z(i+1) = exp(x(i)^3/3);
+    printf('%f    %f    %f\n',x(i),y(i),z(i+1))
+end
+
+printf('\nNote that the exact solution is calculated from the analytical solution y = exp(x^3/3)\n')
+\ No newline at end of file
diff --git a/260/CH13/EX13.3/13_3.sce b/260/CH13/EX13.3/13_3.sce
new file mode 100644
index 000000000..497a84f91
--- /dev/null
+++ b/260/CH13/EX13.3/13_3.sce
@@ -0,0 +1,50 @@
+//Eg-13.3
+//pg-525
+
+clear
+clc
+close()
+x(1) = 0;
+yb(1) = 0;
+y(1) = 1;    //Initial condition
+h = 0.1;
+deff('out = func(in1,in2)','out = in1^2*in2')
+
+//Taking the exact values from the previous problem
+z(1) = exp(x(1)^3/3);
+
+for(i = 2:11)
+    
+    x(i) = 0.1+(i-2)*0.1;
+    z(i+1) = exp(x(i)^3/3);
+end
+
+for(i = 1:10)
+    yb(i+1) = y(i) + h*func(x(i),y(i));
+    y(i+1) = y(i) + h/2*(func(x(i),y(i)) + func(x(i+1),yb(i+1)));
+end
+
+
+
+printf('   x          yPc        yExact\n')
+
+for(i = 2:11)
+    printf('%f    %f    %f\n',x(i),y(i),z(i+1))
+end
+
+for(k = 1:10)
+    a(k) = x(k+1);
+    b(k) = y(k+1);
+    c(k) = z(k+2);
+end
+
+
+clf()
+//plot(a,[b c])
+plot(a,b,'-')
+
+plot(a,c,'.')
+legend('Exact','Euler')
+xlabel('x')
+ylabel('y')
+printf('\n\nA comparison of Eulers method, the predictor-corrector method, and the exact solution\nare presented in the image. As expected, the predictor-corrector method produces a more\naccurate solution than the simple Eulers method.\n')
+\ No newline at end of file
diff --git a/260/CH13/EX13.4/13_4.sce b/260/CH13/EX13.4/13_4.sce
new file mode 100644
index 000000000..f67c1c3a0
--- /dev/null
+++ b/260/CH13/EX13.4/13_4.sce
@@ -0,0 +1,38 @@
+// about the analytical integration for the analytical solution
+//Eg-13.4
+//pg-529
+
+clear
+clc
+
+
+x(1) = 0;
+yb(1) = 0;
+y(1) = 1;    //Initial condition
+h = 0.1;
+deff('out = func(in1,in2)','out = 0.5*(1+in1)*in2^2')
+
+//Taking the exact values using the expression calculated analytically as y = 4/(4-2*x-x^2)
+z(1) = exp(x(1)^3/3);
+
+for(i = 2:11)
+    
+    x(i) = 0.1+(i-2)*0.1;
+    z(i+1) = 4/(4-2*x(i)-x(i)^2);
+end
+
+for(i = 1:10)
+    yb(i+1) = y(i) + h*func(x(i),y(i));
+    y(i+1) = y(i) + h/2*(func(x(i),y(i)) + func(x(i+1),yb(i+1)));
+end
+
+
+
+printf('   x         yPc        yExact\n')
+
+for(i = 2:11)
+    printf('%f    %f    %f\n',x(i),y(i),z(i+1))
+end
+
+printf('\nSolving analytically gives :\n')
+printf('y = 4/(4-2*x-x^2)\n')
+\ No newline at end of file
diff --git a/260/CH13/EX13.5/13_5.sce b/260/CH13/EX13.5/13_5.sce
new file mode 100644
index 000000000..789721088
--- /dev/null
+++ b/260/CH13/EX13.5/13_5.sce
@@ -0,0 +1,34 @@
+//Eg-13.5
+//pg-530
+
+clear
+clc
+
+y(1) = 1;
+x(1) = 0;
+h = 0.1;
+
+deff('out = func(in1,in2)','out = in1^0.5 + in2^0.5')
+
+for(i = 2:11)
+    
+    x(i) = 0.1+(i-2)*0.1;
+    
+end
+
+for(i = 1:10)
+    yb(i+1) = y(i) + h*func(x(i),y(i));
+    y(i+1) = y(i) + h/2*(func(x(i),y(i)) + func(x(i+1),yb(i+1)));
+    bb(i+1) = y(i) + h/2*func(x(i),y(i));
+    b(i+1) = y(i) + h*func(x(i)+h/2,bb(i+1));
+end
+
+
+
+printf('   x         yMP        yPC\n')
+
+for(i = 2:11)
+    printf('%f    %f    %f\n',x(i),b(i),y(i))
+end
+
+printf('The third column presents the results obtained from the method used in Eg-13_4')
+\ No newline at end of file
diff --git a/260/CH13/EX13.6/13_6.sce b/260/CH13/EX13.6/13_6.sce
new file mode 100644
index 000000000..c0fea21a9
--- /dev/null
+++ b/260/CH13/EX13.6/13_6.sce
@@ -0,0 +1,43 @@
+//Eg-13.6
+//pg-533
+
+clear
+clc
+
+//Equation [27] is used for computation.
+
+deff('out = func(in1,in2)','out = in2-in1')
+h = 0.1;
+y(1) = 0;
+
+//The given equation dy/dx = y-x is of the form dy/dx + Py = Q
+// Finally the analytical expression is y = x-exp(x) + 1
+
+//The index starts from 1 and goes up to 11 instead of 0 to 10
+
+for(i = 1:11)
+    x(i) = 0 + (i-1)*h;
+    z(i) = x(i) - exp(x(i)) + 1;
+    
+end
+
+a = (2^0.5-1)/2;
+b = (2-2^0.5)/2;
+c = -(2^0.5)/2;
+d = 1 + (2^0.5)/2;
+
+printf('i      x(i)        y(i)          k1           k2          k3          k4           y(i+1)      y(i+1)Exact\n')
+for(i = 1:10)
+    k1(i) = h*func(x(i),y(i));
+    k2(i) = h*func(x(i)+h/2,y(i)+k1(i)/2);
+    k3(i) = h*func(x(i)+h/2,y(i)+a*k1(i)+b*k2(i));
+    k4(i) = h*func(x(i)+h,y(i)+c*k2(i)+d*k3(i));
+    y(i+1) = y(i) + 1/6*(k1(i)+2*b*k2(i)+2*d*k3(i)+k4(i));
+    
+    printf(' %d    %f    %f    %f    %f    %f    %f    %f    %f\n',i,x(i),y(i),k1(i),k2(i),k3(i),k4(i),y(i+1),z(i+1))
+end
+
+printf(' 11    %f    %f\n',x(i+1),y(i+1))
+printf('\n Therefore, y(1) = %f\n',y(i+1))
+
+printf('\n Refer to the textbook for the analytical solution\n')
+\ No newline at end of file
diff --git a/260/CH13/EX13.7/13_7.sce b/260/CH13/EX13.7/13_7.sce
new file mode 100644
index 000000000..b9b1dec95
--- /dev/null
+++ b/260/CH13/EX13.7/13_7.sce
@@ -0,0 +1,49 @@
+//Eg-13.7
+//pg-537
+
+clear
+clc
+
+A = 0.9;
+B = 0.09;
+y(1) = 1;
+
+deff('out = func(in1,in2)','out = A*in2 - B*in2^2')
+
+h = 0.5;
+//Given the expression of analytical solution : y(t) = 10/(1+9*exp(-0.9*t))
+
+//The index again is 1-11 instead of 0-10
+
+for(i = 1:11)
+    
+    x(i) = 0 + (i-1)*h;
+    yex(i) = 10/(1+9*exp(-0.9*x(i)))
+end
+
+a = (2^0.5-1)/2;
+b = (2-2^0.5)/2;
+c = -(2^0.5)/2;
+d = 1 + (2^0.5)/2;
+
+printf('   x         yRKG       yExact\n')
+for(i = 1:10)
+    k1(i) = h*func(x(i),y(i));
+    k2(i) = h*func(x(i)+h/2,y(i)+k1(i)/2);
+    k3(i) = h*func(x(i)+h/2,y(i)+a*k1(i)+b*k2(i));
+    k4(i) = h*func(x(i)+h,y(i)+c*k2(i)+d*k3(i));
+    y(i+1) = y(i) + 1/6*(k1(i)+2*b*k2(i)+2*d*k3(i)+k4(i));
+    
+    printf('%f    %f    %f\n',x(i+1),y(i+1),yex(i+1))
+end
+
+printf('\nTherefore, it is observed that the RKG solution closely matches \nwith the analytical solution.\n')
+
+
+
+
+
+
+
+
+
diff --git a/260/CH13/EX13.8/13_8.sce b/260/CH13/EX13.8/13_8.sce
new file mode 100644
index 000000000..4d10fae11
--- /dev/null
+++ b/260/CH13/EX13.8/13_8.sce
@@ -0,0 +1,57 @@
+//Eg-13.8
+//pg-542
+
+clear
+clc
+close()
+x(1) = 0;
+x_ending = 200;
+l = 1000;    //no. of divisions
+h = (x_ending-x(1))/l;
+
+y5(1) = 0.01;
+y4(1) = 0.01;
+y(1) = 0.01;
+des = 10^(-5);
+deff('out = F(in1,in2)','out = in2^2*(1-in2)')
+//h = 0.5;
+for(i = 1:150)
+
+k1 = h*F(x(i),y(i));
+k2 = h*F(x(i) + 1/5*h,y(i)+1/5*k1);
+k3 = h*F(x(i) + 3/10*h,y(i)+3/40*k1 + 9/40*k2);
+k4 = h*F(x(i)+3/5*h,y(i) + 3/10*k1 - 9/10*k2 + 6/5*k3);
+k5 = h*F(x(i)+h, y(i - 11/54*k1 + 5/2*k2 - 70/27*k3 + 35/27*k4));
+k6 = h*F(x(i)+7/8*h, y(i) + 1631/55296*k1 + 175/512*k2 + 575/13824*k3 + 44275/110592*k4 + 253/4096*k5);
+
+y5(i+1) = y(i) + (37/378*k1 + 250/621*k3 + 125/594*k4 + 512/1771*k6);
+
+y4(i+1) = y(i) + (2825/27648*k1 + 18575/48384*k3 + 13525/55296*k4 + 277/14336*k5 + 1/4*k6);
+
+err = abs(y5(i+1) - y4(i+1));
+
+if(err > des)
+    n = 0.25;
+elseif(err < des)
+    n = 0.2;
+end
+
+h = h*(abs(des/err))^n;
+
+y(i+1) = y(i) + (37/378*k1 + 250/621*k3 + 125/594*k4 + 512/1771*k6);
+
+x(i+1) = x(i)+h;
+
+end 
+
+printf('The values calculated using RK Fehlberg method are tabulated below\n')
+printf('    x          y\n')
+for(i = 1:150)
+    printf('%f    %f\n',x(i),y(i))
+end
+
+plot(x,y)
+xlabel('t')
+ylabel('y')
+
+//There can be a difference in values this code gives and the code in the textbook gives because the author didn't mention the initial value of h that was used.Final value of x has been taken as 200 and 1000 divisions among will give h = (200-0)/1000.
+\ No newline at end of file
diff --git a/260/CH13/EX13.9/13_9.sce b/260/CH13/EX13.9/13_9.sce
new file mode 100644
index 000000000..33bd15066
--- /dev/null
+++ b/260/CH13/EX13.9/13_9.sce
@@ -0,0 +1,81 @@
+//Eg-13.9
+//pg-546
+
+clear
+clc
+close()
+
+//Note that the subscripts of the variables y have been increased by 1 since the subscript 0 is not possible in scilab!
+
+deff('out = f1(in1,in2,in3)','out = -0.08*in1^0.5 - 2*in1^0.2*in2')
+
+deff('out = f2(in1,in2,in3)','out = -3.5*10^-6*in1^0.2*in2 + 1.6*10^-6*in3^0.3')
+
+deff('out = f3(in1,in2,in3)','out = 2*in1^0.2*in2 - 0.16*in3^0.3')
+
+i = 1;
+k1(1,1) = 0;
+k2(1,1) = 0;
+k3(1,1) = 0;
+k4(1,1) = 0;
+
+
+y1(1) = 0.95;
+y2(1) = 0.05;
+y3(1) = 0;
+
+ti = 0;
+tf = 7;
+l = 1000;
+
+t(1) = 0;
+
+h = (tf-ti)/l;
+
+n = 1 + (tf-ti)/h;
+
+for(i = 2:n)
+    t(i) = t(i-1) + h;
+    k1(i,1) = f1(y1(i-1),y2(i-1),y3(i-1));
+    k1(i,2) = f2(y1(i-1),y2(i-1),y3(i-1));
+    k1(i,3) = f3(y1(i-1),y2(i-1),y3(i-1));
+    
+    y1(i) = y1(i-1) + k1(i,1)*h/2;
+    y2(i) = y2(i-1) + k1(i,2)*h/2;
+    y3(i) = y3(i-1) + k1(i,3)*h/2;
+    k2(i,1) = f1(y1(i),y2(i),y3(i));
+    k2(i,2) = f2(y1(i),y2(i),y3(i));
+    k2(i,3) = f3(y1(i),y2(i),y3(i));
+    
+    y1(i) = y1(i-1) + k2(i,1)*h/2;
+    y2(i) = y2(i-1) + k2(i,2)*h/2;
+    y3(i) = y3(i-1) + k2(i,3)*h/2;
+    k3(i,1) = f1(y1(i),y2(i),y3(i));
+    k3(i,2) = f2(y1(i),y2(i),y3(i));
+    k3(i,3) = f3(y1(i),y2(i),y3(i));
+    
+    y1(i) = y1(i-1) + k3(i,1)*h;
+    y2(i) = y2(i-1) + k3(i,2)*h;
+    y3(i) = y3(i-1) + k3(i,3)*h;
+    k4(i,1) = f1(y1(i),y2(i),y3(i));
+    k4(i,2) = f2(y1(i),y2(i),y3(i));
+    k4(i,3) = f3(y1(i),y2(i),y3(i));
+    
+    y1(i) = y1(i-1) + 1/6*h*(k1(i,1) + 2*(k2(i,1)+k3(i,1)) + k4(i,1));
+    y2(i) = y2(i-1) + 1/6*h*(k1(i,2) + 2*(k2(i,2)+k3(i,2)) + k4(i,2));
+    y3(i) = y3(i-1) + 1/6*h*(k1(i,3) + 2*(k2(i,3)+k3(i,3)) + k4(i,3));
+    
+    
+    
+end
+printf('   t          y1          y2          y3\n')
+for(i = 1:100:n)
+    printf('%f    %f    %f    %f\n',t(i),y1(i),y2(i),y3(i))
+end
+
+plot(t,y1,t,y2,t,y3)
+xlabel('t')
+ylabel('y1,y2,y3')
+legend('y1','y2','y3')
+
+//note that the graph shown on the text book uses y1,y2 from o to 1 and y3 from 0 to 0.1
+\ No newline at end of file