Numerical Methods: Mathematica Notebooks
(c) John H. Mathews, 1996

To accompany the text:

Numerical Methods:  for Mathematics, 
Science, and Engineering, 2nd Ed, 1992
Prentice Hall, Inc.
Englewood Cliffs, New Jersey 07632  U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions:      ISBN 0-13-625047-5

This free software is compliments of the author.
E-mail address:  in%"mathews@fullerton.edu"
http://titan.fullerton.edu/~mathews/

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

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

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 = (b - a)/m //N;
  T = Table[0,{m+1}];
  Y = Table[0,{m+1}];
  T[[1]] = a //N;
  Y[[1]] = y0 //N;
  For[k=1, k<=m,k++,
    Y[[k+1]] = Y[[k]] + h f[T[[k]],Y[[k]]] //N;
    T[[k+1]] = a + h k //N];
  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

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

Compare with Mathematica's procedure NDSolve.

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

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

DSolve[{y'[t]==(t-y[t])/2,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) = y0  over  [a,b]  by computing
yk+1  =  yk  +  h/2  [f(tk,yk) + f(tk+1,yk + h f(tk,yk))]    for k = 0,1,...,m-1.
Section	9.3	Heun's Method   Page  441

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 = (b - a)/m //N;
  T = Table[0,{m+1}];
  Y = Table[0,{m+1}];
  T[[1]] = a //N;
  Y[[1]] = y0 //N;
  For[j=1, j<=m,j++,
    K1 = f[T[[j]],Y[[j]]] //N;
    P = Y[[j]] + h K1 //N;
    T[[j+1]] = a + h j //N;
    K2 = f[T[[j+1]],P] //N;
    Y[[j+1]] = Y[[j]] + h (K1 + K2)/2 //N];
  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

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

Compare with Mathematica's procedure NDSolve.

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

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

DSolve[{y'[t]==(t-y[t])/2,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) = y0  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

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=D[Y[t],{t,2}]; 
  Y3=D[Y[t],{t,3}]; 
  Y4=D[Y[t],{t,4}]; 
  Set@@{d1[t_,y_],f[t,y]};
  Set@@{d2[t_,y_],Together[Y2/.{F->f,Y[t]->y}]};
  Set@@{d3[t_,y_],Together[Y3/.{F->f,Y[t]->y}]};
  Set@@{d4[t_,y_],Together[Y4/.{F->f,Y[t]->y}]};
  h = (b - a)/m //N;
  T = Table[0,{m+1}];
  Y = Table[0,{m+1}];
  T[[1]] = a //N;
  Y[[1]] = y0 //N;
  For[j=1, j<=m,j++,
    tj = T[[j]] //N;
    yj = Y[[j]] //N;
    T[[j+1]] = a  + h j //N;
    Y[[j+1]] = yj + h(d1[tj,yj] + h(d2[tj,yj]/2 +
               h(d3[tj,yj]/6 + h d4[tj,yj]/24))); ];
  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

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

Compare with Mathematica's procedure NDSolve.

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

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

DSolve[{y'[t]==(t-y[t])/2,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) = y0
over [a,b]  by using the formula  yk+1 = yk + h/6 [K1 + 2K2 + 2K3 + K4].
Section	9.5	Runge-Kutta Methods   Page  460

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 = (b - a)/m //N;
  T = Table[0,{m+1}];
  Y = Table[0,{m+1}];
  T[[1]] = a //N;
  Y[[1]] = y0 //N;
  For[j=1, j<=m,j++,
    TJ = T[[j]];
    YJ = Y[[j]];
    K1 = h f[TJ,YJ] //N;
    K2 = h f[TJ + h/2, YJ + K1/2] //N;
    K3 = h f[TJ + h/2, YJ + K2/2] //N;
    K4 = h f[TJ + h, YJ + K3] //N;
    Y[[j+1]] = YJ + (K1 + 2K2 + 2K3 + K4)/6 //N;
    T[[j+1]] = a + h j //N];
  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

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

Compare with Mathematica's procedure NDSolve.

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

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

DSolve[{y'[t]==(t-y[t])/2,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) = y0  over  [a,b]  with an error control and variable step-size method.
Section	9.5	Runge-Kutta Methods   Page  461

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 = (b - a)/m //N;
  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

RKFSet = RKF45[f[t,y],0,1.4,0,7,0.0001];
ListPlot[RKFSet,PlotJoined->True,
         PlotRange->{0,6}];
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}];

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) = y0  over [a,b]  by using the predictor	
   pk+1  =  yk  +  h/24  [- 9fk-3  +  37fk-2  -  59fk-1  +  55fk]
and the corrector	
   yk+1  =  yk  +  h/24  [fk-2  -  5fk-1  +  19fk  +  9fk+1].
Section	9.6	Predictor-Corrector Method   Page  471

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 = (b - a)/m //N;
  T = Table[0,{m+1}];
  Y = Table[0,{m+1}];
  T[[1]] = a //N;
  Y[[1]] = y0 //N;
  For[j=1, j<=3,j++,
    TJ = T[[j]];
    YJ = Y[[j]];
    K1 = h f[TJ,YJ] //N;
    K2 = h f[TJ + h/2, YJ + K1/2] //N;
    K3 = h f[TJ + h/2, YJ + K2/2] //N;
    K4 = h f[TJ + h, YJ + K3] //N;
    Y[[j+1]] = YJ + (K1 + 2K2 + 2K3 + K4)/6 //N;
    T[[j+1]] = a + h j //N];
  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 (-9f0 + 37f1 - 59f2 + 55f3);
    T[[k+1]] = a + h k;
    f4 = f[T[[k+1]],P];
    Y[[k+1]] = YK + h2 (f1 - 5f2 + 19f3 + 9f4);
    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 - 5y  with y(0) = 0  over [0,10].

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

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

Compare with Mathematica's procedure NDSolve.

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

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

DSolve[{y'[t]==30-5y[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) = y0  
over [a,b] by using the predictor
   pk+1  =  yk-3  +  4h/3  [2fk-2  -  fk-1  +  2fk]
and the corrector
   yk+1  =  yk-1  +  h/3  [fk-1  +  4fk  +  fk+1].
Section	9.6	Predictor-Corrector Method   Page  472

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 = (b - a)/m //N;
  T = Table[0,{m+1}];
  Y = Table[0,{m+1}];
  T[[1]] = a //N;
  Y[[1]] = y0 //N;
  For[j=1, j<=3,j++,
    TJ = T[[j]];
    YJ = Y[[j]];
    K1 = h f[TJ,YJ] //N;
    K2 = h f[TJ + h/2, YJ + K1/2] //N;
    K3 = h f[TJ + h/2, YJ + K2/2] //N;
    K4 = h f[TJ + h, YJ + K3] //N;
    Y[[j+1]] = YJ + (K1 + 2K2 + 2K3 + K4)/6 //N;
    T[[j+1]] = a + h j //N];
  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]] + 4h (2f1 - f2 + 2f3)/3;
     Pmod = Pnew + 28(Yold-Pold)/29;
     T[[k+1]] = a + h k;
     f4 = f[T[[k+1]],Pmod];
     Y[[k+1]] = Y[[k-1]] + h (f2 + 4f3 + f4)/3;
     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 - 5y  with y(0) = 0  over [0,10].

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

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

Compare with Mathematica's procedure NDSolve.

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

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

DSolve[{y'[t]==30-5y[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) = y0  
over [a,b] by using the predictor	
   pk+1  =  yk-3  +  4h/3  [2fk-2  -  fk-1  +  2fk]
and the corrector
   yk+1  =  (- yk-2Ê+Ê9yk)/8  + 3h/8  [- fk-1  +  2fk  +  fk+1].
Section	9.6	Predictor-Corrector Method   Page  473

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 = (b - a)/m //N;
  T = Table[0,{m+1}];
  Y = Table[0,{m+1}];
  T[[1]] = a //N;
  Y[[1]] = y0 //N;
  For[j=1, j<=3,j++,
    TJ = T[[j]];
    YJ = Y[[j]];
    K1 = h f[TJ,YJ] //N;
    K2 = h f[TJ + h/2, YJ + K1/2] //N;
    K3 = h f[TJ + h/2, YJ + K2/2] //N;
    K4 = h f[TJ + h, YJ + K3] //N;
    Y[[j+1]] = YJ + (K1 + 2K2 + 2K3 + K4)/6 //N;
    T[[j+1]] = a + h j //N];
  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]] + 4h (2f1 - f2 + 2f3)/3;
    Pmod = Pnew + 112(Cold - Pold)/121;
    T[[k+1]] = a + h k;
    f4 = f[T[[k+1]],Pmod];
    Cnew = (9Y[[k]] - Y[[k-2]] + 3h (-f2 + 2f3 + f4))/8;
    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 - 5y  with y(0) = 0  over [0,10].

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

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

Compare with Mathematica's procedure NDSolve.

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

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

DSolve[{y'[t]==30-5y[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) = a and  x(b) = 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) = 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)  +  [(b-u(b))/v(b)] v(t).

First, execute the following group of cells:

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 = (b-a)/m//N;
    T = Table[0,{m+1}];
    X = Table[0,{m+1}];
    Y = Table[0,{m+1}];
    T[[0+1]] = a//N;
    X[[0+1]] = alpha//N;
    Y[[0+1]] = beta//N;
    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 + (k1+2k2+2k3+k4)/6;
        Y[[j+1]] = yj + (r1+2r2+2r3+r4)/6;
        T[[j+1]] = a + h*j//N; ]];
    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[[m+1,2]] sol2[[i,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/(4t^2)) x= t^(1/2) 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/(4t^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}},
         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 t^(1/2));
g[t_] = f[t]/.{C->-2.960284381/4,C1->0.2354853843/4}

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"}];

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) = a  and  x(b) = b
over [a,b] by using the finite difference method of order O(h2). 
Section	9.9	Finite-Difference Method   Page  496

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

First, execute the following group of cells:

Off[General::spell1];
Clear[CoeffMat,a,b,alpha,beta,n];
CoeffMat[a_,b_,alpha_,beta_,n_] :=
Module[{h},
  h  = (b-a)/(n)//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 + h/2 p[Vt[[1]]])  alpha;
  Vb[[n-1]] = Vb[[n-1]] + (1 - h/2 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 - h/2 p[Vt[[j+1]]]; ]];
  Clear[j];
  For[j=1,j<=n-2,j++,
    Module[{},
    Vc[[j]] = -1 + h/2 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/(4t^2)) x= t^(1/2) 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/(4t^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}},
         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 t^(1/2));
g[t_] = f[t]/.{C->-2.960284381/4,C1->0.2354853843/4}

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"}];