Numerical Methods: Mathematica 3.0 Notebooks
(c) John H. Mathews, 1998
To accompany the text:
Numerical Methods:  for Mathematics, 
Science, and Engineering, 2nd Ed, 1992
Prentice Hall, Inc.
Simon & Schuster,  One Lake Street
Upper Saddle River, New Jersey 07458  U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions:      ISBN 0-13-625047-5
www.Prenhall.com
This free software is compliments of the author.
E-mail address:      in%"mathews@fullerton.edu"

CONTENTS

Chapter 9.   Solution of Differential Equations
	Algorithm 9.1    Euler's Method.
	Algorithm 9.2    Heun's Method.
	Algorithm 9.3    Taylor's Method of Order 4.
	Algorithm 9.4    Runge-Kutta Method of Order 4.
	Algorithm 9.5    Runge-Kutta-Fehlberg Method RKF45.
	Algorithm 9.6    Adams-Bashforth-Moulton Method.
	Algorithm 9.7    Milne-Simpson Method.
	Algorithm 9.8    The Hamming Method.
	Algorithm 9.9    Linear Shooting Method.
	Algorithm 9.10  Finite-Difference Method.

Clear["Global`*"];
Off[General::"spell1"];
Clear[a, alpha, a0, ABM, b, beta, b0, c0, CoeffMat, d0, d1, d2, d3, 
        d4, Euler, f0, f1, f2, f3, F0, FindDDE, h, Hamming, Heun, i, j, k, 
        K1, K2, K3, K4, LinearBdd, m, Milne, n, n0, p, p0, q, q0, r, r0, 
        RKbdd4, RKF45, Runge, tol, T, TaylorMeth, TJ, TriMat, va, vd, vc, vb, 
        x, X, y0, Y, YJ];
On[General::"spell1"];

Algorithm 9.1  Euler's Method.  To approximate the solution of the initial value problem  y' =f (t, y)  with   y (a) = (Subscript[y, 0])  over  [a, b]  by computing  (Subscript[y, k + 1])  =  (Subscript[y, k]) + h f ((Subscript[t, k]), (Subscript[y, k]))  for  k = 0, 1, ⋯, m - 1.
Section  9.2   Euler's Method   Page  435

Execute this cell to activate Euler's Method.
Euler[f0_,a_,b_,y0_,m_]
f0_ is the function f(t,y),
a_ is the left endpoint,
b_ is the right endpoint,
y0_is the initial value y(a),
m_ is the number of subintervals.

Off[General::"spell1"]; Clear[Euler, f0, a, b, y0, m, T, Y]; \(Euler[f0_, a_, b_, y0_, m_] := 
	
      Module[{f, h, k}, 
		Set@@{f[t_, y_], f0}; 
		
        h = N[\(b - a\)\/m]; 
		T = Table[0, {m + 1}]; 
		
        Y = Table[0, {m + 1}]; 
		
        T[[1]] = N[a]; 
		
        Y[[1]] = N[y0]; 
		
        For[k = 1, k ≤ m, \(k++\), 
          Y[[k + 1]] = 
            N[Y[[k]] + 
                h\ f[T[[k]], 
                    Y[[k]]]]; 
			
          T[[k + 1]] = N[a + h\ k]]; 

        		sol = Transpose[{T, Y}]; 
		Print[" "]; 
		
        Print["The D.E. is  y' = ", f[t, x]]; 
		
        Print["Euler's method is used with ", m, " subintervals."]; 

        		Return[sol]]; 
On[General::"spell1"]; \)

Chapter 9, Example 9.4, Page 433.  Solve  y' =((t - y)/2)  with  y (0) = 1  over  [0, 3].

Clear[f, t, y]; f[t_, y_] = ((t - y)/2.0);
Print["
y' = ",f[t, y]];

EulerSet = Euler[f[t, y], 0, 3, 1, 12]; \(ListPlot[EulerSet, PlotJoined -> True, PlotRange -> {0, 2}, 
      PlotStyle -> RGBColor[1, 0, 0]]; 
sol\)

Compare with Mathematica's procedure NDSolve.

NDSolve[{y'[t]==(1/2) (t - y[t]),y[0] == 1},y[t],{t, 0, 3}];Plot[Evaluate[y[t] /. %], {t, 0, 3}, PlotRange -> {0, 2}, 
      PlotStyle -> RGBColor[0, 1, 0]];

Or use Mathematica's procedure DSolve to get an analytic solution, in this case.

DSolve[{y'[t]==(1/2) (t - y[t]),y[0] == 1},y[t],t]

Algorithm 9.2  Heun's Method.  To approximate the solution of the initial value problem  y' =f (t, y)  with   y (a) = (Subscript[y, 0])  over  [a, b]  by computing  (Subscript[y, k + 1]) = (Subscript[y, k])+((h )/2)[f ((Subscript[t, k]), (Subscript[y, k])) + 
              f ((Subscript[t, k + 1]), (Subscript[y, k]) + h f ((Subscript[t, k]), (Subscript[y, k])))]  for  k = 0, 1, ⋯, m - 1.
Section  9.3   Heun's Method   Page  441

Execute this cell to activate Heun's Method.
Heun[f0_,a_,b_,y0_,m_]
f0_ is the function f(t,y),
a_ is the left endpoint,
b_ is the right endpoint,
y0_is the initial value y(a),
m_ is the number of subintervals.

Off[General::"spell1"]; Clear[Heun, a, b, y0, m, T, Y]; \(Heun[f0_, a_, b_, y0_, m_] := 
	
      Module[{f, h, j}, 
		Set@@{f[t_, y_], f0}; 
		
        h = N[\(b - a\)\/m]; 
		T = Table[0, {m + 1}]; 
		
        Y = Table[0, {m + 1}]; 
		
        T[[1]] = N[a]; 
		
        Y[[1]] = N[y0]; 
		
        For[j = 1, j ≤ m, \(j++\), 
			
          K1 = N[f[T[[j]], 
                Y[[j]]]]; 
			
          P = N[Y[[j]] + h\ K1]; 
			
          T[[j + 1]] = N[a + h\ j]; 

          			K2 = 
            N[f[T[[j + 1]], P]]; 
			
          Y[[j + 1]] = 
            N[Y[[j]] + 
                \(h\/2\) \((K1 + K2)\)]]; 
		sol = Transpose[{T, Y}]; 

        		Print[" "]; 
		Print["The D.E. is  y' = ", f[t, x]]; 

        		Print["Heun's method is used with ", m, 
          " subintervals."]; 
		Return[sol]]; 

    \(On[General::"spell1"]; \)\)

Chapter 9, Example 9.6, Page 439.  Solve  y' =((t - y)/2)  with  y (0) = 1  over  [0, 3].

Clear[f, t, y]; f[t_, y_] = ((t - y)/2.0);
Print["
y' = ",f[t, y]];

HeunSet = Heun[f[t, y], 0, 3, 1, 12]; \(ListPlot[HeunSet, PlotJoined -> True, PlotRange -> {0, 2}, 
      PlotStyle -> RGBColor[1, 0, 0]]; 
sol\)

Compare with Mathematica's procedure NDSolve.

NDSolve[{y'[t]==(1/2) (t - y[t]),y[0] == 1},y[t],{t, 0, 3}];Plot[Evaluate[y[t] /. %], {t, 0, 3}, PlotRange -> {0, 2}, 
      PlotStyle -> RGBColor[0, 1, 0]];

Or use Mathematica's procedure DSolve to get an analytic solution, in this case.

DSolve[{y'[t]==(1/2) (t - y[t]),y[0] == 1},y[t],t]

Algorithm 9.3  Taylor's Method of Order 4.  To approximate the solution of the initial value problem  y' =f (t, y)  with   y (a) = (Subscript[y, 0])  over  [a, b]  by evaluating  (y')'  ,  ((y')')'   and  ((y')')'′   and using the Taylor polynomial at each step.
Section  9.4   Taylor Series Method   Page  448

Execute this cell to activate Taylor's Method of Order 4.
TaylorMeth[f0_,a_,b_,y0_,m_]
f0_ is the function f(t,y),
a_ is the left endpoint,
b_ is the right endpoint,
y0_is the initial value y(a),
m_ is the number of subintervals.

\(Off[General::"spell1"]; 
\(Clear[TaylorMeth, f0, a, b, y0, m, T, Y]; 
    \)\)TaylorMeth[f0_, a_, b_, y0_, m_]:=
	Module[{d1, d2, d3, d4, f, h, j},
		Set@@{f[t_, y_], f0};
		Y'[t]=F[t, Y[t]];
		\(Y2 = ∂\_{t, 2}Y[t]\);
		\(Y3 = ∂\_{t, 3}Y[t]\);
		\(Y4 = ∂\_{t, 4}Y[t]\);
		Set@@{d1[t_, y_], f[t, y]};
		Set@@{d2[t_, y_], 
                  Together[ReplaceAll[Y2, {F -> f, Y[t] -> y}]]};
		Set@@{d3[t_, y_], 
                  Together[ReplaceAll[Y3, {F -> f, Y[t] -> y}]]};
		Set@@{d4[t_, y_], 
                  Together[ReplaceAll[Y4, {F -> f, Y[t] -> y}]]};
		h = N[((b - a)/m)];
		T = Table[0, {m + 1}];
		Y = Table[0, {m + 1}];
		T[[1]] = N[a];
		Y[[1]] = N[y0];
		\(For[j = 1, j ≤ m, \(j++\), 
			
                tj = N[T[[j]]]; 
			
                yj = N[Y[[j]]]; 
			
                T[[j + 1]] = 
                  N[a + h\ j]; 
			
                Y[[j + 1]] = 
                  yj + h\ \((
                        d1[tj, yj] + 
                          h\ \((1\/2\ d2[tj, yj] + 
                                h\ \((1\/6\ d3[tj, yj] + 1\/24\ h\ d4[tj, yj])
                                    \))\))\); ]\);
		sol = Transpose[{T, Y}];
		Print[" "];
		Print["The D.E. is  y' = ", f[t, x]];
		Print["Taylor's method of order 4 is used with ", m, 
                " subintervals."];
		Return[sol]];
On[General::"spell1"];

Chapter 9, Example 9.8, Page 445.  Solve  y' =((t - y)/2)  with  y (0) = 1  over  [0, 3].

Clear[f, t, y]; f[t_, y_] = ((t - y)/2.0);
Print["
y' = ",f[t, y]];

TaylorSet = TaylorMeth[f[t, y], 0, 3, 1, 12]; \(ListPlot[TaylorSet, PlotJoined -> True, PlotRange -> {0, 2}, 
      PlotStyle -> RGBColor[1, 0, 0]]; 
sol\)

Compare with Mathematica's procedure NDSolve.

NDSolve[{y'[t]==(1/2) (t - y[t]),y[0] == 1},y[t],{t, 0, 3}];Plot[Evaluate[y[t] /. %], {t, 0, 3}, PlotRange -> {0, 2}, 
      PlotStyle -> RGBColor[0, 1, 0]];

Or use Mathematica's procedure DSolve to get an analytic solution, in this case.

DSolve[{y'[t]==(1/2) (t - y[t]),y[0] == 1},y[t],t]

Algorithm 9.4 Runge-Kutta Method of Order 4.  To approximate the solution of the initial value problem  y' =f (t, y)  with   y (a) = (Subscript[y, 0])  over  [a, b]   by using the formula  (Subscript[y, k + 1]) = (Subscript[y, k])+((h )/6)[(Subscript[K, 1]) + 2  (Subscript[K, 2]) + 2  (Subscript[K, 3]) + (Subscript[K, 4])].
Section  9.5   Runge-Kutta Methods   Page  460

Execute this cell to activate Runge-Kutta Method of Order 4.
Runge[f0_,a_,b_,y0_,m_]
f0_ is the function f(t,y),
a_ is the left endpoint,
b_ is the right endpoint,
y0_is the initial value y(a),
m_ is the number of subintervals.

Off[General::"spell1"]; Clear[Runge, f0, a, b, y0, m, T, Y]; \(Runge[f0_, a_, b_, y0_, m_] := 
	
      Module[{f, h, j, K1, K2, K3, K4, TJ, YJ}, 
		Set@@{f[t_, y_], f0}; 

        		h = N[\(b - a\)\/m]; 
		T = Table[0, {m + 1}]; 
		
        Y = Table[0, {m + 1}]; 
		
        T[[1]] = N[a]; 
		
        Y[[1]] = N[y0]; 
		
        For[j = 1, j ≤ m, \(j++\), 
			
          TJ = T[[j]]; 
			
          YJ = Y[[j]]; 
			
          K1 = N[h\ f[TJ, YJ]]; 
			K2 = N[h\ f[TJ + h\/2, YJ + K1\/2]]; 
          
			K3 = N[h\ f[TJ + h\/2, YJ + K2\/2]]; 
			
          K4 = N[h\ f[TJ + h, YJ + K3]]; 
			
          Y[[j + 1]] = 
            N[YJ + 1\/6\ \((K1 + 2\ K2 + 2\ K3 + K4)\)]; 
			
          T[[j + 1]] = N[a + h\ j]]; 

        		Print[" "]; 
		Print["The D.E. is  y' = ", f[t, x]]; 

        		Print["The Runge-Kutta method is used with ", m, 
          " subintervals."]; 
		sol = Transpose[{T, Y}]; 
		
        Return[sol]]; 
\(On[General::"spell1"]; \)\)

Chapter 9, Example 9.10, Page 454.  Solve  y' =((t - y)/2)  with  y(0) = 1  over  [0, 3].

Clear[f, t, y];
f[t_, y_] = ((t - y)/2.0);
Print["
y' = ",f[t, y]];

RungeSet = Runge[f[t, y], 0, 3, 1, 12]; \(ListPlot[RungeSet, PlotJoined -> True, PlotRange -> {0, 2}, 
      PlotStyle -> RGBColor[1, 0, 0]]; 
sol\)

Compare with Mathematica's procedure NDSolve.

NDSolve[{y'[t]==(1/2) (t - y[t]),y[0] == 1},y[t],{t, 0, 3}];Plot[Evaluate[y[t] /. %], {t, 0, 3}, PlotRange -> {0, 2}, 
      PlotStyle -> RGBColor[0, 1, 0]];

Or use Mathematica's procedure DSolve to get an analytic solution, in this case.

DSolve[{y'[t]==(1/2) (t - y[t]),y[0] == 1},y[t],t]

Algorithm 9.5  Runge-Kutta-Fehlberg Method RKF45.  To approximate the solution of the initial value problem  y' =f (t, y)  with   y (a) = (Subscript[y, 0])  over  [a, b]  with an error control and variable step-size method.
Section  9.5   Runge-Kutta Methods   Page  461

Execute this cell to activate Runge-Kutta-Fehlberg Method RKF45.
RKF45[f0_,a_,b_,y0_,m_,tol_]
f0_ is the function f(t,y),
a_ is the left endpoint,
b_ is the right endpoint,
y0_is the initial value y(a),
m_ is the number of subintervals,
tol_ is the tolerance.

Off[General::"spell1"]; Clear[RKF45, f0, a, b, y0, m, tol, T, Y]; \(RKF45[f0_, a_, b_, y0_, m_, tol_] := 
	
      Module[{f, h, j}, 
		Set@@{f[t_, y_], f0}; 
		a2 = 1\/4; 
		
        b2 = 1\/4; 
		a3 = 3\/8; 
		b3 = 3\/32; 
		c3 = 9\/32; 
		
        a4 = 12\/13; 
		b4 = 1932\/2197; 
		c4 = \(-\(7200\/2197\)\); 

        		d4 = 7296\/2197; 
		a5 = 1; 
		b5 = 439\/216; 
		
        c5 = \(-8\); 
		d5 = 3680\/513; 
		e5 = \(-\(845\/4104\)\); 

        		a6 = 1\/2; 
		b6 = \(-\(8\/27\)\); 
		c6 = 2; 
		
        d6 = \(-\(3544\/2565\)\); 
		e6 = 1859\/4104; 
		
        f6 = \(-\(11\/40\)\); 
		r1 = 1\/360; 
		
        r3 = \(-\(128\/4275\)\); 
		r4 = \(-\(2197\/75240\)\); 
		
        r5 = 1\/50; 
		r6 = 2\/55; 
		n1 = 25\/216; 
		
        n3 = 1408\/2565; 
		n4 = 2197\/4104; 
		n5 = \(-\(1\/5\)\); 

        		h = N[\(b - a\)\/m]; 
		T = {}; 
		T = AppendTo[T, N[a]]; 

        		Y = {}; 
		Y = AppendTo[Y, N[y0]]; 
		hmin = h\/64; 
		
        hmax = 64\ h; 
		max1 = 200; 
		j = 1; 
		
        tj = T[[1]]; 
		
        br = b - 0.00001\ Abs[b]; 
		
        While[T[[j]] < b, 
          If[T[[j]] + h > br, 
            h = b - T[[j]]]; 
			
          tj = T[[j]]; 
			
          yj = Y[[j]]; 
			y1 = yj; 
          
			k1 = h\ f[tj, y1]; 
			y2 = yj + b2\ k1; 
			
          k2 = h\ f[tj + a2\ h, y2]; 
			y3 = yj + b3\ k1 + c3\ k2; 

          			k3 = h\ f[tj + a3\ h, y3]; 
			
          y4 = yj + b4\ k1 + c4\ k2 + d4\ k3; 
			
          k4 = h\ f[tj + a4\ h, y4]; 
			
          y5 = yj + b5\ k1 + c5\ k2 + d5\ k3 + e5\ k4; 
			
          k5 = h\ f[tj + a5\ h, y5]; 
			
          y6 = yj + b6\ k1 + c6\ k2 + d6\ k3 + e6\ k4 + f6\ k5; 
			
          k6 = h\ f[tj + a6\ h, y6]; 
			
          err = Abs[r1\ k1 + r3\ k3 + r4\ k4 + r5\ k5 + r6\ k6]; 
			
          ynew = yj + n1\ k1 + n3\ k3 + n4\ k4 + n5\ k5; 
			
          If[err < tol || h < 2\ hmin, 
				
            \(Module[{}, 
					Y = AppendTo[Y, ynew]; 
					
              If[tj + h > br, 
						T = AppendTo[T, N[b]], 

                						T = AppendTo[T, N[tj + h]]]; 
					
              j = j + 1; 
					
              tj = T[[j]]; ]; \)]; 

          			If[err == 0, s = 0, s = 0.84\ \((\(tol\ h\)\/err)\)\^0.25]; 

          			If[s < 0.75 && h > 2\ hmin, h = h\/2]; 
			
          If[s > 1.50 && 2\ h < hmax, h = 2\ h]; 
			
          If[max1 == j, Break[]]; ]; 
		sol = Transpose[{T, Y}]; 
		
        Print[" "]; 
		Print["The D.E. is  y' = ", f[t, x]]; 
		
        Print["The Runge-Kutta-Fehlberg method is used."]; 
		
        Return[sol]]; 
\(On[General::"spell1"]; \)\)

Chapter 9, Example 9.12, Page 459.  Solve  y' = 1 + (y^2)  with  y (0) = 0  over  [0, 1.4].

Clear[f, t, y]; f[t_, y_] = 1 + (y^2);
Print["
y' = ",f[t, y]];

RKFSet = RKF45[f[t, y], 0, 1.4, 0, 7, 0.0001]; \(ListPlot[RKFSet, PlotJoined -> True, PlotRange -> {0, 6}, 
      PlotStyle -> RGBColor[1, 0, 0]]; 
sol\)

Compare with Mathematica's procedure NDSolve.

NDSolve[{y'[t]==1 + y[t]^2,y[0] == 0},y[t],{t, 0, 1.4}];Plot[Evaluate[y[t] /. %], {t, 0, 1.4}, PlotRange -> {0, 6}, 
      PlotStyle -> RGBColor[0, 1, 0]];

Or use Mathematica's procedure DSolve to get an analytic solution, in this case.

DSolve[{y'[t]==1 + y[t]^2,y[0] == 0},y[t],t]

Algorithm 9.6 Adams-Bashforth-Moulton Method.  To approximate the solution of the initial value problem  y' =f (t, y)  with   y (a) = (Subscript[y, 0])  over  [a, b]   by using the predictor  
	(Subscript[p, k + 1]) = (Subscript[y, k])+((h )/24)[- 9 (Subscript[f, k - 3]) + 37  (Subscript[f, k - 2]) - 59  (Subscript[f, k - 1]) + 
              55  (Subscript[f, k])]
and the corrector
	(Subscript[y, k + 1]) = (Subscript[y, k])+((h )/24)[(Subscript[f, k - 2]) - 5  (Subscript[f, k - 1]) + 19  (Subscript[f, k]) + 9  (Subscript[f, k + 1])].
Section  9.6   Predictor-Corrector Method   Page  471

Execute this cell to activate Adams-Bashforth-Moulton Method.
ABM[F0_,a_,b_,y0_,m_]
f0_ is the function f(t,y),
a_ is the left endpoint,
b_ is the right endpoint,
y0_is the initial value y(a),
m_ is the number of subintervals.

Off[General::"spell1"]; Clear[ABM, F0, a, b, y0, m]; \(ABM[F0_, a_, b_, y0_, m_] := 
	
      Module[{f, f0, f1, f2, f3, h, h2, j, k}, 
		Set@@{f[t_, y_], F0}; 

        		h = N[\(b - a\)\/m]; 
		T = Table[0, {m + 1}]; 
		
        Y = Table[0, {m + 1}]; 
		
        T[[1]] = N[a]; 
		
        Y[[1]] = N[y0]; 
		
        For[j = 1, j ≤ 3, \(j++\), 
			
          TJ = T[[j]]; 
			
          YJ = Y[[j]]; 
			
          K1 = N[h\ f[TJ, YJ]]; 
			K2 = N[h\ f[TJ + h\/2, YJ + K1\/2]]; 
          
			K3 = N[h\ f[TJ + h\/2, YJ + K2\/2]]; 
			
          K4 = N[h\ f[TJ + h, YJ + K3]]; 
			
          Y[[j + 1]] = 
            N[YJ + 1\/6\ \((K1 + 2\ K2 + 2\ K3 + K4)\)]; 
			
          T[[j + 1]] = N[a + h\ j]]; 

        		f0 = 
          f[T[[1]], 
            Y[[1]]]; 
		
        f1 = f[T[[2]], 
            Y[[2]]]; 
		
        f2 = f[T[[3]], 
            Y[[3]]]; 
		
        f3 = f[T[[4]], 
            Y[[4]]]; 
		h2 = h\/24; 

        		For[k = 4, k ≤ m, \(k++\), 
			
          TK = T[[k]]; 
			
          YK = Y[[k]]; 
			
          P = YK + h2\ \((\(-9\)\ f0 + 37\ f1 - 59\ f2 + 55\ f3)\); 
			
          T[[k + 1]] = a + h\ k; 
			
          f4 = f[T[[k + 1]], P]; 
			
          Y[[k + 1]] = 
            YK + h2\ \((f1 - 5\ f2 + 19\ f3 + 9\ f4)\); 
			f0 = f1; 

          			f1 = f2; 
			f2 = f3; 
			
          f3 = f[T[[k + 1]], 
              Y[[k + 1]]]]; 
		
        sol = Transpose[{T, Y}]; 
		Print[" "]; 
		
        Print["The D.E. is  y' = ", f[t, x]]; 
		
        Print["The Adams-Bashforth-Moulton method is used with ", m, 
          " subintervals."]; 
		Return[sol]]; 

    \(On[General::"spell1"]; \)\)

Chapter 9, Example 9.14, Page 470.  Solve  y' = 30 - 5  y  with  y (0) = 0  over  [0, 10].

Clear[f, t, y]; f[t_, y_] = 30 - 5 y;
Print["
y' = ",f[t, y]];

ABMSet = ABM[f[t, y], 0, 10, 1, 65]; ListPlot[ABMSet, PlotJoined -> True, 
      PlotRange -> {{0, 10.2}, {0, 7.4}}, 
      PlotStyle -> RGBColor[1, 0, 0], 
      Ticks -> {Range[0, 10, 1], Range[0, 7, 1]}, 
      AxesLabel -> {"t", "y"}];

ABMSet = ABM[f[t, y], 0, 10, 1, 37]; ListPlot[ABMSet, PlotJoined -> True, 
      PlotRange -> {{0, 10.2}, {0, 7.4}}, 
      PlotStyle -> RGBColor[0, 0, 1], 
      Ticks -> {Range[0, 10, 1], Range[0, 7, 1]}, 
      AxesLabel -> {"t", "y"}];

Compare with Mathematica's procedure NDSolve.

NDSolve[{y'[t]==30 - 5 y[t],y[0] == 0},y[t],{t, 0, 10}];Plot[Evaluate[y[t] /. %], {t, 0, 10}, 
      PlotRange -> {{0, 10.2}, {0, 7.4}}, 
      PlotStyle -> RGBColor[0, 1, 0]];

Or use Mathematica's procedure DSolve to get an analytic solution, in this case.

DSolve[{y'[t]==30 - 5 y[t],y[0] == 0},y[t],t]

Algorithm 9.7 Milne-Simpson Method. To approximate the solution of the initial value problem  y' =f (t, y)  with   y (a) = (Subscript[y, 0])  over  [a, b]   by using the predictor
	(Subscript[p, k + 1]) = (Subscript[y, k - 3])+((4 h )/3)[2  (Subscript[f, k - 2]) - (Subscript[f, k - 1]) + 2  (Subscript[f, k])]
and the corrector
	(Subscript[y, k + 1]) = (Subscript[y, k - 1])+((h )/3)[(Subscript[f, k - 1]) + 4  (Subscript[f, k]) + (Subscript[f, k + 1])].
Section  9.6   Predictor-Corrector Method   Page  472

Execute this cell to activate Milne-Simpson Method.
Milne[F0_,a_,b_,y0_,m_]
f0_ is the function f(t,y),
a_ is the left endpoint,
b_ is the right endpoint,
y0_is the initial value y(a),
m_ is the number of subintervals.

Off[General::"spell1"]; Clear[Milne, F0, a, b, y0, m]; \(Milne[F0_, a_, b_, y0_, m_] := 
	
      Module[{f, f0, f1, f2, f3, h, j, k}, 
		Set@@{f[t_, y_], F0}; 
		
        h = N[\(b - a\)\/m]; 
		T = Table[0, {m + 1}]; 
		
        Y = Table[0, {m + 1}]; 
		
        T[[1]] = N[a]; 
		
        Y[[1]] = N[y0]; 
		
        For[j = 1, j ≤ 3, \(j++\), 
			
          TJ = T[[j]]; 
			
          YJ = Y[[j]]; 
			
          K1 = N[h\ f[TJ, YJ]]; 
			K2 = N[h\ f[TJ + h\/2, YJ + K1\/2]]; 
          
			K3 = N[h\ f[TJ + h\/2, YJ + K2\/2]]; 
			
          K4 = N[h\ f[TJ + h, YJ + K3]]; 
			
          Y[[j + 1]] = 
            N[YJ + 1\/6\ \((K1 + 2\ K2 + 2\ K3 + K4)\)]; 
			
          T[[j + 1]] = N[a + h\ j]]; 

        		f1 = 
          f[T[[2]], 
            Y[[2]]]; 
		
        f2 = f[T[[3]], 
            Y[[3]]]; 
		
        f3 = f[T[[4]], 
            Y[[4]]]; 
		Pold = 0; 

        		Yold = 0; 
		
        For[k = 4, k ≤ m, \(k++\), 
			
          Pnew = Y[[k - 3]] + 
              4\/3\ h\ \((2\ f1 - f2 + 2\ f3)\); 
			
          Pmod = Pnew + \(28\ \((Yold - Pold)\)\)\/29; 
			
          T[[k + 1]] = a + h\ k; 
			
          f4 = f[T[[k + 1]], Pmod]; 

          			Y[[k + 1]] = 
            Y[[k - 1]] + 
              1\/3\ h\ \((f2 + 4\ f3 + f4)\); 
			Pold = Pnew; 
			
          Yold = Y[[k + 1]]; 
			
          f1 = f2; 
			f2 = f3; 
			
          f3 = f[T[[k + 1]], 
              Y[[k + 1]]]]; 
		
        sol = Transpose[{T, Y}]; 
		Print[" "]; 
		
        Print["The D.E. is  y' = ", f[t, x]]; 
		
        Print["The Milne-Simpson method is used with ", m, 
          " subintervals."]; 
		Return[sol]]; 

    \(On[General::"spell1"]; \)\)

Chapter 9, Example 9.14, Page 470.  Solve  y' = 30 - 5  y  with  y (0) = 0  over  [0, 10].

Clear[f, t, y]; f[t_, y_] = 30 - 5 y;
Print["
y' = ",f[t, y]];

MilneSet = Milne[f[t, y], 0, 10, 1, 110]; ListPlot[MilneSet, PlotJoined -> True, 
      PlotRange -> {{0, 10.2}, {0, 7.4}}, 
      PlotStyle -> RGBColor[1, 0, 0], 
      Ticks -> {Range[0, 10, 1], Range[0, 7, 1]}, 
      AxesLabel -> {"t", "y"}];

MilneSet = Milne[f[t, y], 0, 10, 1, 93]; ListPlot[MilneSet, PlotJoined -> True, 
      PlotRange -> {{0, 10.2}, {0, 7.4}}, 
      PlotStyle -> RGBColor[0, 0, 1], 
      Ticks -> {Range[0, 10, 1], Range[0, 7, 1]}, 
      AxesLabel -> {"t", "y"}];

Compare with Mathematica's procedure NDSolve.

NDSolve[{y'[t]==30 - 5 y[t],y[0] == 0},y[t],{t, 0, 10}];Plot[Evaluate[y[t] /. %], {t, 0, 10}, 
      PlotRange -> {{0, 10.2}, {0, 7.4}}, 
      PlotStyle -> RGBColor[0, 1, 0]];

Or use Mathematica's procedure DSolve to get an analytic solution,in this case.

DSolve[{y'[t]==30 - 5 y[t],y[0] == 0},y[t],t]

Algorithm  9.8  The Hamming Method.  To approximate the solution of the initial value problem  y' =f (t, y)  with   y (a) = (Subscript[y, 0])  over  [a, b]   by using the predictor
	(Subscript[p, k + 1]) = (Subscript[y, k - 3])+((4 h )/3)[2  (Subscript[f, k - 2]) - (Subscript[f, k - 1]) + 2  (Subscript[f, k])]	
and the corrector
	(Subscript[y, k + 1]) =(( -(Subscript[y, k - 2]) + 9 (Subscript[y, k]))/8)+((3 h )/8)[-(Subscript[f, k - 1]) + 2  (Subscript[f, k]) + (Subscript[f, k + 1])].
Section  9.6   Predictor-Corrector Method   Page  473

Execute this cell to activate The Hamming Method.
Hamming[F0_,a_,b_,y0_,m_]
f0_ is the function f(t,y),
a_ is the left endpoint,
b_ is the right endpoint,
y0_is the initial value y(a),
m_ is the number of subintervals.

Off[General::"spell1"]; Clear[Hamming, F0, a, b, y0, m]; \(Hamming[F0_, a_, b_, y0_, m_] := 
	
      Module[{f, f0, f1, f2, f3, h, j, k}, 
		Set@@{f[t_, y_], F0}; 
		
        h = N[\(b - a\)\/m]; 
		T = Table[0, {m + 1}]; 
		
        Y = Table[0, {m + 1}]; 
		
        T[[1]] = N[a]; 
		
        Y[[1]] = N[y0]; 
		
        For[j = 1, j ≤ 3, \(j++\), 
			
          TJ = T[[j]]; 
			
          YJ = Y[[j]]; 
			
          K1 = N[h\ f[TJ, YJ]]; 
			K2 = N[h\ f[TJ + h\/2, YJ + K1\/2]]; 
          
			K3 = N[h\ f[TJ + h\/2, YJ + K2\/2]]; 
			
          K4 = N[h\ f[TJ + h, YJ + K3]]; 
			
          Y[[j + 1]] = 
            N[YJ + 1\/6\ \((K1 + 2\ K2 + 2\ K3 + K4)\)]; 
			
          T[[j + 1]] = N[a + h\ j]]; 

        		f1 = 
          f[T[[2]], 
            Y[[2]]]; 
		
        f2 = f[T[[3]], 
            Y[[3]]]; 
		
        f3 = f[T[[4]], 
            Y[[4]]]; 
		Pold = 0; 

        		Cold = 0; 
		
        For[k = 4, k ≤ m, \(k++\), 
			
          TK = T[[k]]; 
			
          YK = Y[[k]]; 
			
          Pnew = Y[[k - 3]] + 
              4\/3\ h\ \((2\ f1 - f2 + 2\ f3)\); 
			
          Pmod = Pnew + \(112\ \((Cold - Pold)\)\)\/121; 
			
          T[[k + 1]] = a + h\ k; 
			
          f4 = f[T[[k + 1]], Pmod]; 

          			Cnew = 
            1\/8\ \((
                9\ Y[[k]] - 
                  Y[[k - 2]] + 
                  3\ h\ \((\(-f2\) + 2\ f3 + f4)\))\); 
			
          Y[[k + 1]] = 
            Cnew + \(9\ \((Pnew - Cnew)\)\)\/121; 
			Pold = Pnew; 

          			Cold = Cnew; 
			f1 = f2; 
			f2 = f3; 
			
          f3 = f[T[[k + 1]], 
              Y[[k + 1]]]]; 
		
        sol = Transpose[{T, Y}]; 
		Print[" "]; 
		
        Print["The D.E. is  y' = ", f[t, x]]; 
		
        Print["Hamming's method is used with ", m, " subintervals."]; 
        
		Return[sol]]; 
\(On[General::"spell1"]; \)\)

Chapter 9, Example 9.14, Page 470.  Solve  y' = 30 - 5  y  with  y (0) = 0  over  [0, 10].

Clear[f, t, y]; f[t_, y_] = 30 - 5 y;
Print["
y' = ",f[t, y]];

HammingSet = Hamming[f[t, y], 0, 10, 1, 70]; ListPlot[HammingSet, PlotJoined -> True, 
      PlotRange -> {{0, 10.2}, {0, 7.4}}, 
      PlotStyle -> RGBColor[1, 0, 0], 
      Ticks -> {Range[0, 10, 1], Range[0, 7, 1]}, 
      AxesLabel -> {"t", "y"}];

HammingSet = Hamming[f[t, y], 0, 10, 1, 50]; ListPlot[HammingSet, PlotJoined -> True, 
      PlotRange -> {{0, 10.2}, {0, 7.4}}, 
      PlotStyle -> RGBColor[0, 0, 1], 
      Ticks -> {Range[0, 10, 1], Range[0, 7, 1]}, 
      AxesLabel -> {"t", "y"}];

Compare with Mathematica's procedure NDSolve.

NDSolve[{y'[t]==30 - 5 y[t],y[0] == 0},y[t],{t, 0, 10}];Plot[Evaluate[y[t] /. %], {t, 0, 10}, 
      PlotRange -> {{0, 10.2}, {0, 7.4}}, 
      PlotStyle -> RGBColor[0, 1, 0]];

Or use Mathematica's procedure DSolve to get an analytic solution, in this case.

DSolve[{y'[t]==30 - 5 y[t],y[0] == 0},y[t],t]

Algorithm 9.9  Linear Shooting Method.   To approximate the solution of the boundary value problem   (x' )' = p(t) x' (t)+q (t) x (t)+r (t)  with  x (a) = α  and  x (b) = β  over the interval  [a, b]   by the Runge-Kutta method of order 4.
Section   9.8   Boundary Value Problems   Page  488

The method involves solving a two systems of equations over  [a, b]
u' = y (t)						with	u (a) = α
y' = p (t) y (t) + q (t) u (t) + r (t)	and	y (a)=u'(a)=0.
then
v' = z (t)						with	v (a) = 0 
z' = p (t) z (t) + q (t) v (t) + r (t)	and	z (a)=v'(a)=1.
 Finally, the desired solution x(t) is the linear combination
 x (t) = u (t)+((β - u (b))/(v (b)))v(t).

Execute this cell to activate Linear Shooting Method.
RKbdd4[f0_,a_,b_,alpha_,beta_,m_]
LinearBdd[p0_,q0_,r0_,a_,b_,alpha_,beta_,m_]
f0_ is the function f(t,y),
p0_ is the function p(t),
q0_ is the function q(t),
r0_ is the function r(t),
a_ is the left endpoint,
b_ is the right endpoint,
alpha_is the initial value y(a),
beta_is the initial value y(b),
m_ is the number of subintervals.

Off[General::"spell1"]; Clear[RKbdd4, f0, a, b, alpha, beta, m, T, X, Y]; \(\(RKbdd4[f0_, a_, b_, alpha_, beta_, m_] := 
	
      Module[{f, h, i, j}, 
		Set@@{f[t_, x_, y_], f0}; 
		
        h = N[\(b - a\)\/m]; 
		T = Table[0, {m + 1}]; 
		
        X = Table[0, {m + 1}]; 
		Y = Table[0, {m + 1}]; 
		
        T[[0 + 1]] = N[a]; 
		
        X[[0 + 1]] = N[alpha]; 
		
        Y[[0 + 1]] = N[beta]; 
		
        For[j = 1, j ≤ m, \(j++\), 
			
          Module[{}, 
				
            tj = T[[j]]; 
				
            xj = X[[j]]; 
				
            yj = Y[[j]]; 
				
            k1 = h\ yj; 
				r1 = h\ f[tj, xj, yj]; 
				
            k2 = h\ \((yj + r1\/2)\); 
				
            r2 = h\ f[tj + h\/2, xj + k1\/2, yj + r1\/2]; 
				
            k3 = h\ \((yj + r2\/2)\); 
				
            r3 = h\ f[tj + h\/2, xj + k2\/2, yj + r2\/2]; 
				
            k4 = h\ \((yj + r3)\); 
				
            r4 = h\ f[tj + h, xj + k3, yj + r3]; 
				
            X[[j + 1]] = 
              xj + 1\/6\ \((k1 + 2\ k2 + 2\ k3 + k4)\); 
				
            Y[[j + 1]] = 
              yj + 1\/6\ \((r1 + 2\ r2 + 2\ r3 + r4)\); 
				
            T[[j + 1]] = 
              N[a + h\ j]; ]]; 
		sol0 = Transpose[{T, X}]; 
		
        Return[sol0]]; \)\)

Clear[LinearBdd, p0, q0, r0, a, b, alpha, beta, m]; LinearBdd[p0_, q0_, r0_, a_, b_, alpha_, beta_, m_]:=
	Module[{f1, f2, p, q, r},
		Set@@{p[t_], p0};
		Set@@{q[t_], q0};
		Set@@{r[t_], r0};
		Set@@{f1[t_, x_, y_], p[t] y + q[t] x + r[t]};
		Set@@{f2[t_, x_, y_], p[t] y + q[t] x};
		sol1 = RKbdd4[f1[t, x, y], a, b, alpha, 0, m];
		sol2 = RKbdd4[f2[t, x, y], a, b, 0, 1, m];
		Clear[i];
		sol = Table[{
                    sol1[[i, 1]], 
                    sol1[[i, 2]] + 
                      (((beta - 
                              sol1[[m + 1, 2
                                ]]) 
                          sol2[[i, 2
                            ]])/sol2)
                          [[m + 1, 2
                          ]]}, {i, 1, m + 1}];
		Print[" "];
		Print["The D.E. is  x'' = ",f1[t, x, y]/.y->x'];
		Print["The Linear Shooting method is used with ", 
                m, " subintervals."];
		Return[sol];];
On[General::"spell1"];

Chapter 9, Exercise 7, Page 489.  Solve   (x' )' +(1/t) x' +(1-(1/(4 (t^2))))x = (Sqrt[t]) cos (t)  over  [1, 6]  with  x (1) = 1.0  and  x (6) = -0.5.

Clear[p, q, r, t]; p[t_] = -(1/t); q[t_] = -(1 - (1/(4 (t^2)))); r[t_] = (Sqrt[t]) Cos[t];

BddValSet = LinearBdd[p[t], q[t], r[t], 1, 6, 1.0, -0.5, 25]; 
    \(ListPlot[BddValSet, PlotJoined -> True, 
      PlotRange -> {{0, 6.2}, {\(-3.2\), 1.2}}, 
      PlotStyle -> RGBColor[1, 0, 0], 
      Ticks -> {Range[0, 6, 1], Range[\(-3\), 2, 1]}, 
      AxesLabel -> {"t", "x"}]; 
sol\)

Compare with the analytic solution to the D.E.

f[t_] = 
      ((t Cos[t] + (t^2) Sin[t] - 4 C Sin[t] + 4 C1 Cos[t])/(4 
          (Sqrt[t]))); \(g[t_] = 
      ReplaceAll[
        f[t], {C -> \(-\(2.960284381\/4\)\), 
          C1 -> 0.2354853843\/4}]; 
Print["g[t] = ", g[t]]; \)

\(exactsol = 
      N[Table[{1.0 + \(i - 1\)\/5, g[1.0 + \(i - 1\)\/5]}, {i, 1, 26}], 8]; 

    ListPlot[exactsol, PlotJoined -> True, 
      PlotRange -> {{0, 6.2}, {\(-3.2\), 1.2}}, 
      Ticks -> {Range[0, 6, 1], Range[\(-3\), 2, 1]}, 
      AxesLabel -> {"t", "x"}, 
      PlotStyle -> RGBColor[0, 1, 0]]; \)

Algorithm 9.10.  Finite-Difference Method.  To approximate the solution of the boundary value problem   (x' )' = p(t) x' (t)+ q (t) x (t)+r (t)  with  x (a) = α  and  x (b) = β  over the interval  [a, b]  by using the finite difference method of order  O((h^2)). 
Section  9.9   Finite-Difference Method   Page  496

Set up the linear system to be solved. It is formed by the subroutine CoeffMat[n].

Execute this cell to activate Finite-Difference Method.
CoeffMat[a_,b_,alpha_,beta_,n_]
TriMat[va_,vd_,vc_,vb_,n_]
FineDDE[a_,b_,alpha_,beta_,n_]
a_ is the left endpoint,
b_ is the right endpoint,
alpha_is the initial value y(a),
beta_is the initial value y(b),
n_ is the number of subintervals.

Off[General::"spell1"]; Clear[CoeffMat, a, b, alpha, beta, n]; \(\(CoeffMat[a_, b_, alpha_, beta_, n_] := 
	
      Module[{h}, 
		h = N[\(b - a\)\/n]; 
		Clear[j]; 
		
        Vt = Table[a + h\ j, {j, n - 1}]; 
		Va = Table[0, {n - 2}]; 
		
        Vc = Table[0, {n - 2}]; 
		Vd = Table[1, {n - 1}]; 
		
        Vb = Table[0, {n - 1}]; 
		Clear[j]; 
		
        For[j = 1, j ≤ n - 1, \(j++\), 
			
          Module[{}, 
            \(Vb[[j]] = 
              \(-h\^2\)\ r[Vt[[j]]]; 
            \)]]; 
		
        Vb[[1]] = 
          Vb[[1]] + 
            \((1 + 1\/2\ h\ p[Vt[[1]]])
                \)\ alpha; 
		
        Vb[[n - 1]] = 
          Vb[[n - 1]] + 
            \((1 - 1\/2\ h\ 
                    p[Vt[[n - 1]]])\)\ 
              beta; 
		Clear[j]; 
		
        For[j = 1, j ≤ n - 1, \(j++\), 
			
          Module[{}, 
            \(Vd[[j]] = 
              2 + h\^2\ q[Vt[[j]]]; 
            \)]]; 
		Clear[j]; 
		
        For[j = 1, j ≤ n - 2, \(j++\), 
			
          Module[{}, 
            \(Va[[j]] = 
              \(-1\) - 
                1\/2\ h\ p[Vt[[j + 1]]]; 
            \)]]; 
		Clear[j]; 
		
        For[j = 1, j ≤ n - 2, \(j++\), 
			
          Module[{}, 
            \(Vc[[j]] = 
              \(-1\) + 
                1\/2\ h\ p[Vt[[j]]]; 
            \)]]; ]; \)\)

Clear[TriMat, va, vd, vc, vb, n]; \(\(TriMat[va_, vd_, vc_, vb_, n_] := 
	
      Module[{a0, b0, c0, d0, n0, k, tem}, 
		n0 = n - 1; 
		a0 = va; 

        		b0 = vb; 
		c0 = vc; 
		d0 = vd; 
		Clear[x0]; 
		
        x0 = Table[0, {n0}]; 
		
        For[k = 2, k ≤ n0, \(k++\), 
			
          Module[{}, 
				
            tem = a0[[k - 1]]\/d0
                  [[k - 1]]; 
				
            d0[[k]] = 
              d0[[k]] - 
                tem\ c0[[k - 1]]; 

            				b0[[k]] = 
              b0[[k]] - 
                tem\ b0[[k - 1]]]]; 

        		x0[[n0]] = 
          b0[[n0]]\/d0
              [[n0]]; 
		
        For[k = n0 - 1, 1 ≤ k, \(k--\), 
			
          x0[[k]] = 
            \(b0[[k]] - 
                c0[[k]]\ 
                  x0[[k + 1]]\)\/d0
                [[k]]]; 
		
        Return[x0]]; \)\)

Clear[FindDDE, a, b, alpha, beta, n]; FineDDE[a_, b_, alpha_, beta_, n_]:=
	Module[{i},
		CoeffMat[a, b, alpha, beta, n];
		TriMat[Va, Vd, Vc, Vb, n];
		sol = Table[{0.0, 0.0}, {n + 1}];
		Clear[i];
		\(For[i = 1, i ≤ n - 1, \(i++\), 
			
                \(sol[[i + 1]] = {
                    Vt[[i]], 
                    x0[[i]]}; \)]\);
		sol[[1]] = {a, alpha};
		sol[[n + 1]] = {b, beta}
                  ;
		Print[" "];
		Print["The D.E. is  x'' = ",p[t] y + q[t] x + r[t]/.y->x'];
		Print[
                "The Finite-Difference method is used with ", m, 
                " subintervals."];
		Return[sol];];
On[General::"spell1"];

Chapter 9, Exercise 7, Page 497.  Solve   (x' )' +(1/t) x' +(1-(1/(4 (t^2))))x = (Sqrt[t]) cos (t)  over  [1, 6]  with  x (1) = 1.0  and  x (6) = -0.5.

Clear[p, q, r, t]; p[t_] = -(1/t); q[t_] = -(1 - (1/(4 (t^2)))); r[t_] = (Sqrt[t]) Cos[t];

FineDESet = FineDDE[1, 6, 1.0, -0.5, 25]; \(ListPlot[FineDESet, PlotJoined -> True, 
      PlotRange -> {{0, 6.2}, {\(-3.2\), 1.2}}, 
      PlotStyle -> RGBColor[1, 0, 0], 
      Ticks -> {Range[0, 6, 1], Range[\(-3\), 2, 1]}, 
      AxesLabel -> {"t", "x"}]; 
sol\)

Compare with the analytic solution to the D.E.

f[t_] = 
      ((t Cos[t] + (t^2) Sin[t] - 4 C Sin[t] + 4 C1 Cos[t])/(4 
          (Sqrt[t]))); \(g[t_] = 
      ReplaceAll[
        f[t], {C -> \(-\(2.960284381\/4\)\), 
          C1 -> 0.2354853843\/4}]; 
Print["g[t] = ", g[t]]; \)

\(exactsol = 
      N[Table[{1.0 + \(i - 1\)\/5, g[1.0 + \(i - 1\)\/5]}, {i, 1, 26}], 8]; 

    ListPlot[exactsol, PlotJoined -> True, 
      PlotRange -> {{0, 6.2}, {\(-3.2\), 1.2}}, 
      Ticks -> {Range[0, 6, 1], Range[\(-3\), 2, 1]}, 
      AxesLabel -> {"t", "x"}, 
      PlotStyle -> RGBColor[0, 1, 0]]; \)