NASA Technical Memorandum 101609 

Solutions of Differential 

Equations with Regular 

Coefficients by The Methods of 

Richmond and Runge-Kutta 

C. R. Cockrell 

August 1989 

- 

(MaSA-lVl-lO 1609) SOLUTIOES OF DIPPEPEITIAL H89-27060 - EQUATIOHS UITB REGU LA8 COEPPICXEITS B? THB 

Langley Research Center) 20 p CSCL 201s Uaclas 

tlmEODS OF RICHMOND AYD BUR6B-KUTTa (UASA. 

63/32 0224 104 

National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23665-5225 
SOLUTIONS OF DIFFERENTIAL EQUATIONS WITH REGULAR COEFFICIENTS 

BY THE METHODS OF RICHMOND AND RUNGE-KUTTA 

ABSTRACT 

Numerical solutions of the differential equation which 

describes the electric field within an inhomogeneous layer of 

permittivity, upon which a perpendicularly-polarized plane wave is 

incident, are considered. Richmond's method and the Runge-Kutta 

method are compared for linear and exponential profiles of 

permittivities. These two approximate solutions are also compared 

with the exact solutions. 
INTRODUCTION 

The treatment of perpendicularly-polarized plane wave 

incidence on inhomogeneous dielectric or plasma layers has been 

investigated by many authors [l-31. Their efforts have been 

directed toward solving the differential equations which describe 

the electromagnetic behavior of fields within these layers. 

Solutions for continuous inhomogeneous layers require, in general, 

numerical techniques. Exact solutions for linear and exponential 

profiles are available [4]. 

Many methods have been used to solve the differential 

equations describing the fields within an inhomogeneous layer. 

Dividing the layer into a number of smaller segments and 

approximating the inhomogeneous profile by constant permittivity 

values within each segment, the multilayer boundary value problem 

is then solved. Approximate integral solutions are also available 

as are WKB asymptotic solutions where the inhomogeneous profile 

(or permittivity) gradient is small [5,4]. For these solutions 

the computation time is too long or the model is too crude to 

yield accurate results. 

Practical solutions have been obtained using the Runge-Kutta 

method [6] and a step-by-step numerical method described by 

Richmond [7]. Both solutions have been programmed on a personal 

computer and the results compared. A comparison with exact 

solutions, linear and exponential, has also been made. 

The primary purpose of this paper is to explain and 

demonstrate the Richmond and Runge-Kutta methods for solving a 

differential equation in which its parameters may or may not be 

known analytically. For definiteness, therefore, the differential 

equation which represents the electric field within a plane 

inhomogeneous layer of dielectric or plasma, upon which a 

perpendicularly-polarized plane wave has been assumed incident, 

has been considered. An efficient solution of this differential 

equation for this ideal model is needed since the same 

differential equation exists for realistic problems such as 

aperture antennas which are bounded by inhomogeneous dielectric or 

plasma layers. 

2 
INHOMOGENEOUS PLANE LAYER-PERPENDICULAR POLARIZATION 

A plane wave is incident on a plane inhomogeneous layer of 

dielectric or plasma as shown in figure 1. The incident electric 

field which is polarized in the x-direction propagates in free 

space (region I) at an angle 8 with the normal to the layer. The 

transmitted electric field which is also x-polarized exists into 

free space (region 111) at an angle 8. The partial differential 

equation for the electric field in the free space regions and in 

the inhomogeneous layer is written as 

2 2 

(1) a Ex(YIZ) a Ex(Y,Z) 2 + 2 + k Cr(Z) Ex(Y,Z) = 0 

aY2 az 

where k is the free space wavenumber and cr(z) is the relative 

dielectric constant of the medium. Since the relative dielectric 

constant in free space (regions I and 111) equals one, the 

solutions of equation (1) in these regions are readily found: 

~I(y,z) = E~[ e jkzcose + 

e -j kzcose] ej kysine 

I11 j kzcose ej kysine Ex (YlZ) = Ete (3) 

where Ei is the intensity of the incident electric field, R is the 

reflection coefficient at the interface, and Et is the transmitted 

electric field intensity. However, in region 11, since E (z) can 

be an arbitrary function of z, the solution of equation (1) must 

be determined numerically, except for some specialized profiles 

'such as linear and exponential. 

Before attempting the solution of equation (1) in region 11, 

it is reduced to an ordinary differential equation. This is 

readily accomplished by the method of separation of variables and 

the continuity of tangential fields throughout the media. The 

continuity requirement forces the y-variation in the fields to be 

the same as given in equations (2) and (3). The differential 

equation describing the z-variation of the electric field is 

written as 

r 

2 I1 

(4) 

2 I1 d Ex (2) 

+ k2( cr(z) - sin e )Ex (z) = 0 

dz2 

3 
Richmond's Method 

One method of solving equation (4) for region I1 is reported in a paper by Richmond [7]. In his method a step-by-step 

numerical procedure is used to find the electric field within the 

layer at successive points. To begin the solution, the field at 

z=O and z=hl which is a small distance within the layer, must be 

determined first. The field at z=O is assumed to be unity and the 

field at z=h is found from the Maclaurin series expansion at 

z=O; i.e., 

3 

h+- 21 Ex(Z)1h3+ . . (5) 3 ! dz3 

d2Ex(Z) I z=o 

I1 Ex (h) = Ex(o)+ dEx(Z) dz /h + 1 d22 

z=o z=o 

where 

Ex(O ) = 1 

= jk cose dz 

z=o 

2 

dz 

3 

2 sin e ) Ex(z)I = - k 2 ( cC- 2 

z=o 

dzEc 1 Ex(z)I 2 = - k2[ jk ( cC- sin 8 ) cose + - 

z=o dZ3 

and cC is the complex relative permittivity at the point z= 0, and 

dz 

point z=O+. Substituting equation (6) into equation (5) yields 

- dcC is the derivative of the permittivity with respect to z at the 

Ex I1 (h) = 1 - 7 kLhL (&+-sine) 2 -z k2h3 E+' 

where E+ is the real relative permittivity, E+' is the derivative 

of the real relative permittivity with respect to z, and tang+ is 

the loss tangent. The subscript + denotes evaluation at z=O+. 

The Taylor series expansion of the electric field about the 

point z=nh is written as 

4 
E1'(z) = Ex I1 (nh) + dE'l(z)l dz ( z - nh) + - 1 d2E'1(z) I ( z -nh) 2 + ... (8) 

dz2 z=nh 

X 

z=nh 

Evaluating equation (8) at z=nh+h and z=nh-h and adding the 

result, thus giving 

E:' (nh+h) I1 

r 2Ex (nh) + h2 

2 2 I1 = k ( cr(z) - sin 8 ) Ex (z) I 

z=nh z=nh dz2 

(9) 

4 I1 

Ex 4 (') I = k4(c$z)-sin2e )E:'(z) 

z=nh dz 

11' I1 - 2k cr' (z) Ex (z) - k2c;' (z) Ex (z) I (11) 

z=nh 

with 

(nh) - E:'(nh-h) I1 

1- 

11' EX 

E.. (nh) = A 11 

Employing the notation of Richmond, 

(10) substituted, the electric field 

2 [ 2 - k 22 h (cn - sin2e) ] En+l 

where ET denotes equation (11). The 

equation (9), with equation 

in the layer becomes 

2 % 

En En-l +h 12 - 

dielectric or plasma layer is 

divided into N segment so that Nh equals D, the thickness of the 

layer. The criteria for choosing the appropriate step size h is 

discussed in reference 7. The solution of equation (13) is 

determined in a step-by-step procedure as n varies from unity to 

N; i.e., for n=l the field at z=2h is computed from equation (13) 

since the fields at z=O and z=h are given by equations (6) and (7) 

and the fourth derivative of the field at z=h is given by 

equations (10)-(12). By proceeding through the layer in this 

manner, the field at the region I and region I1 interface is 

determined. The field at this interface then enables one to 

determine the reflection coefficient. 

5 
Runae-Kutta Method 

The equations needed for solving second-order differential 

equations of the form 

by the Runge-Kutta method are given in reference 6. 

Equation (4) is written in this form as 

L II 

d Ex (z) = - k 2 ( Cr(Z) - sin 2 8 )Ex I1 (z) 

Using the notation adapted earlier, the electric field at a 

general point z=nh where h is the spacing within the layer is 

represented as En and at z=nh+h as En+l. From the fourth-order 

Runge-Kutta method, the solution takes the form 

=En+h 

I I 2 1 

- lk + - k + - k En+l -En+z 1 3 2 6 3 

where the primes denote the derivative with respect to z and 

kl=hf nh E [ n) 

1 h 'h nh + h, En+ En + g kl] 

'h nh + h, En+ h En + 

k2) 

6 
with 

f(nh, En] = - k2[ s(nh) - sin 2 8 )En 

f(nh+ 1 h'h h, En+ - 2 E + - 8 '1 = 

cinh+ z) h - sin28] [ En+ h 'h 

En + kl] 

nh+ h, En+ h En + 

' 

k2) = 

cr(nh+ h) - sin28] [En+ h E + 'h 

k2 I J n 

In this procedure, n varies from zero to N-1. The initial values 

Eo and Eo are given by the first two equations in equation (6). 

I 

7 
DISCUSSION RESULTS 

Linear and exponential profiles were chosen to demonstrate 

the Richmond and Runge-Kutta methods for solving the differential 

equations. A comparison of these methods is made. These two 

profiles were selected since they also have exact solutions as 

shown in the Appendix. 

The particular linear profile chosen is shown in figure 2. 

The relative permittivity cr begins at a value of minus one (the 

exit point) and rises to a value of two at the front interface of 

the inhomogeneous layer. The real and imaginary parts of the 

electric fields at front interface for the Richmond (equation 

(13)) and theo Runge-Kutta (equation (16)) methods are shown in 

Table I for 5 increments in the angle of incidence. Also shown 

are the complex values of the electric fields for the exact 

solutions (Airy-equation (A-8)). The approximate computational 

time for each value in the Richmond method was 8 seconds compared 

to 12 seconds for the Runge-Kutta method. Both methods compare 

very favorably with the exact solutions as can be readily seen in 

Table I. 

TABLE I 

ELECTRIC FIELD DISTRIBUTION AT INTERFACE OF INHOMOGENEOUS 

LAYER-LINEAR PROFILE 

I m= 3k and b=-ll 

Angle 

0 

5 

10 

15 

20 

25 

30 

35 

40 

45 

50 

55 

60 

65 

70 

75 

80 

85 

90 

Richmond 

Real 

0.9918 

0.9954 

1.0063 

1.0243 

1.0488 

1.0794 

1.1150 

1.1550 

1.1983 

1.2435 

1.2894 

1.3345 

1.3775 

1.4166 

1.4508 

1.4788 

1.4995 

1.5124 

1.5167 

Imag 

0.9178 

0.9155 

0.9086 

0.8969 

0.8800 

0.8577 

0.8296 

0.7953 

0.7543 

0.7065 

0.6515 

0.5895 

0.5206 

0.4452 

0.3640 

0.2777 

0.1874 

0.0944 

0.0000 

Runge-Kutta 

Real 

0.9919 

0.9955 

1.0065 

1.0244 

1.0489 

1.0794 

1.1152 

1.1552 

1.1985 

1.2437 

1.2895 

1.3336 

1.3775 

1.4167 

1.4509 

1.4789 

1.4997 

1.5125 

1.5168 

Imag 

0.9179 

0.9156 

0.9086 

0.8969 

0.8800 

0.8577 

0.8296 

0.7953 

0.7543 

0.7065 

0.6516 

0.5895 

0.5206 

0.4452 

0.3640 

0.2777 

0.1874 

0.0944 

0.0000 

Exact-Airy 

Real 

0.9919 

0.9955 

1.0065 

1.0244 

1.0489 

1.0794 

1.1152 

1.1552 

1.1985 

1.2437 

1.2895 

1.3346 

1.3775 

1.4167 

1.4509 

1.4789 

1.4997 

1.5125 

1.5168 

Imag 

0.9179 

0.9156 

0.9086 

0.8969 

0.8800 

0.8577 

0.8296 

0.7953 

0.7543 

0.7065 

0.6516 

0.5895 

0.5206 

0.4452 

0.3640 

0.2777 

0.1874 

0.0944 

0.0000 . 

8 
The exponential profile selected for comparison is shown in 

figure 3. The relative permittivity cr begins at a value of one 

and rises to a value of e at the front interface of the 

inhomogeneous layer. The values of the real and imaginary parts 

of the electric fields at the front interface foF the Richmond and 

Runge-Kutta methods are shown in Table I1 for 5 increments in the 

angle of incidence. Also shown are the cornplFx valuFs of the 

electric fields for two angles of incidence ( 0 and 30 ) for the 

exact solution (Bessel-equation (A-17 )): these two angles were 

selected because, for the parameters used in these computations, 

they yield integer order Bessel functions which are easily 

calculable. The approximate computational time for each value in 

the Richmond method was 9 seconds compared to 12 seconds for the 

Runge-Kutta method. Both methods compare very favorably with each 

other as well as to the two exact solutions as can be readily seen 

in Table 11. 

TABLE I1 

ELECTRIC FIELD DISTRIBUTION AT INTERFACE OF INHOMOGENEOUS 

LAYER-EXPONENTIAL PROFILE 

9ngle 

0 

5 

10 

15 

20 

25 

30 

35 

40 

45 

50 

55 

60 

65 

70 

75 

80 

85 

90 

Richmond 

Real 

0.3736 

0.3766 

0.3851 

0.3993 

0.4187 

0.4427 

0.4711 

0.5028 

0.5372 

0.5731 

0.6096 

0.6456 

0.6798 

0.7112 

0.7385 

0.7610 

0.7776 

0.7878 

0.7912 

Imag 

0.7411 

0.7393 

0.7340 

0.7249 

0.7119 

0.6945 

0.6725 

0.6455 

0.6130 

0.5748 

0.5308 

0.4809 

0.4251 

0.3639 

0.2976 

0.2273 

0.1535 

0.0774 

0.0000 

Runge-Kutta 

Real 

0.3738 

0.3767 

0.3853 

0.3994 

0.4188 

0.4430 

0.4713 

0.5030 

0.5373 

0.5733 

0.6098 

0.6457 

0.6800 

0.7113 

0.7387 

0.7611 

0.7777 

0.7880 

0.7914 

Imag 

0.7411 

0.7393 

0.7340 

0.7249 

0.7119 

0.6945 

0.6725 

0.6455 

0.6130 

0.5749 

0.5308 

0.4809 

0.4252 

0.3640 

0.2978 

0.2273 

0.1535 

0.0774 

0.0000 

Exact-Bessel 

Real 

0.3738 

0.4722 

Imag 

0.7411 

0.6722 

Reflection coefficients at the front interface of the 

inhomogeneous layers are also of interest. The equations required 

for computing the reflection coefficients are given in reference 

7. Magnitude and phase of the reflection coefficients for the 

linear and exponential profiles discussed in this paper are shown 

in Tables I11 and IV. 

9 
TABLE I11 

REFLECTION COEFFICIENTS AT INTERFACE OF INHOMOGENEOUS 

LAYER-LINEAR PROFILE 

I Richmond 

I1 11 

I m= 3k and b=-1 ] 

Runge-Kutta 

Phase 

in radians 

-. 0548 

-.0480 

-.0277 

0.0030 

0.0542 

0.1159 

0.1919 

0.2828 

0.3897 

0.5127 

0.6540 

0.8154 

1.0007 

1.2139 

1.4642 

1.7644 

2.1351 

2.5981 

3.1416 

0 

5 

10 

15 

20 

25 

30 

35 

40 

45 

50 

55 

60 

65 

70 

75 

80 

85 

90 

Magnitude 

0.2929 

0.2936 

0.2957 

0.2992 

0.3043 

0.3112 

0.3201 

0.3314 

0.3457 

0.3638 

0.3870 

0.4169 

0.4564 

0.5092 

0.5811 

0.6788 

0.8048 

0.9367 

1.0000 

0.2925 

0.2932 

0.2954 

0.2991 

0.3044 

0.3116 

0.3208 

0.3324 

0.3473 

0.3659 

0.3896 

0.4202 

0.4604 

0.5137 

0.5862 

0.6840 

0.8090 

0.9384 

1.0000 

Phase 

in radians 

-.0825 

-.0757 

-. 0553 

-.0211 

0.0268 

0.0889 

0.1655 

0.2571 

0.3644 

0.4885 

0.6308 

0.7934 

0.9797 

1.1946 

1.4463 

1.7485 

2.1221 

2.5901 

3.1416 

10 
TABLE IV 

90 11 1.0000 

REFLECTION COEFFICIENTS AT INTERFACE OF INHOMOGENEOUS 

LAYER-EXPONENTIAL PROFILE /I 

-3.1416 

Richmond II 

Magnitude 

0.2825 

0.2837 

0.2874 

0.2935 

0.3024 

0.3144 

0.3297 

0.3489 

0.3727 

0.4021 

0.4382 

0.4828 

0.5376 

0.6051 

0.6969 

0.7823 

0.8828 

0.9663 

1.0000 

Runge-Kutta 

Phase 

in radians 

-2.4801 

-2.4791 

-2.4762 

-2.4714 

-2.4651 

-2.4575 

-2.4494 

-2.4414 

-2.4345 

-2.4302 

-2.4303 

-2.4374 

-2.4553 

-2.4889 

-2.5455 

-2.6338 

-2.7631 

-2.9367 

-3.1416 

Angle 

0 

5 

10 

15 

20 

25 

30 

35 

40 

45 

50 

55 

60 

65 

70 

75 

80 

85 

I Magnitude 

0.2788 

0.2800 

0.2836 

0.2898 

0.2986 

0.3105 

0.3258 

0.3449 

0.3686 

0.3979 

0.4340 

0.4784 

0.5334 

0.6009 

0.6831 

0.7791 

0.8808 

0.9656 

Phase 

in radians 

-2.4876 

-2.4863 

-2.4780 

-2.4710 

-2.4631 

-2.4540 

-2.4450 

-2.4372 

-2.4320 

-2.4311 

-2.4373 

-2.4542 

-2.4868 

-2.5427 

-2.6306 

-2.7601 

-2.9348 

-2.4833 

. 

11 
CONCLUDING REMARKS 

The numerical techniques of Richmond and Runge-Kutta were 

shown to be quite efficient for solving the second-order 

differential equation which describes the field within a plane 

inhomogeneous layer of permittivity, upon which a perpendicularlypolarized plane wave has been assumed incident. Richmond's method 

seems to be little more efficient than the Runge-Kutta method, at 

least for the two relative permittivity profiles considered. 

Although the techniques of Richmond and Runge-Kutta were 

applied to analytical profiles, they should be equally valid for 

arbitrary relative permittivity profiles or even profiles which 

are known discretely. This is important because realistic 

relative permittivity profiles may take on almost any shape. 

12 
APPENDIX 

Exact solutions of equation (4) are available for only a few 

relative permittivity profiles [3,4]. A brief derivation of the 

solutions of equation (4) for linear and exponential profiles is 

given here. 

For the linear profile, the relative permittivity is assumed 

to vary as 

cr(z) = mz + b (A-1) 

where m is the slope of the profile and b is the beginning value 

of the profile ( at z=O ). Substituting equation (A-1) into 

equation (4) yields 

2 

dz2 

d Ex(z) 2 2 + k ( mz + b - sin 8 ) Ex(z) = 0 (A-2 1 

By making 

equation 

d2E"(V) 

the substitution 

z=Mv+B 

A-2) becomes 

(A-3 1 

22 + { M k [ m M v + m B + b - sin28 1) Ex(v) = 0 (A-4 1 A 

dv2 

By equating the term in braces to -v and solving for M and B, it 

is found that 

(A-5) 

where 

(A-6 1 

2 sin 8 - b 

m 

B= 

From equations (A-3) and (A-6), it follows that 

13 
2 

2 k7 

v = - ( mz + b- sin 8 )[ (A-7 1 

The solutions of equation (A-5) are the two independent Airy 

functions, Ai(v) and Bi(v) [4]. The general solution of equation 

(A-5), therefore is written as 

Ex(v) = CIAi(v) + CZBi(v) (A-8 1 

where C1 and C2 are constants. These constants are determined by 

applying the appropriate boundary conditions at z=O. The boundary 

conditions are given by the first two equations in equation (6) 

but which must be expressed in terms of the new independent v. 

The evaluation of the constants is given in matrix form as 

I 

Bi(v0) - Bi(v0) 1 

I (A-9) 1 

- Ai(v0) Ai (VO) 

where 

I I 1 ID = Ai (vo) Bi (vo) - Ai (vo) Bi (vo) 

(A-10) 

and the primes denote the derivative with respect to v. For the 

solutions in this paper m=3k and b=-1. Once the C's are determined, 

the electric field at any point in the layer is computed by using 

equation (8). 

For the exponential profile, the relative permittivity is 

assumed to vary as 

cr(z) = e sz (A-11) 

Substituting equation (A-11) into equation (4) yields 

2 

dz2 

+ k2 [ esz - sin Ex(z) = 0 

d Ex(z) 

Following the substitution given in reference 4, 

(A-12) 

14 
(A-13) 

the differential equation is now written as 

2 

dv2 

2 d Ex(v) 

+v- dEx (v) + ( v2 - p 2 ) Ex(v) = 0 dv V 

where 

If p is zero or a positive integer, the general 

equation (A-14) becomes 

Ex(v) = B J (v) + B Y (v) 1P 2P 

where J (v) is the Bessel function of the first kind, P 

(A-14) 

(A-15) 

solution of 

(A-16) 

Y (v) is the 

P 

Bessel function of the second kind, and the B's are constants. 

Upon applying the appropriate boundary conditions at z=O, the 

matrix for the B's becomes 

where 

I I \ 

ID = J (vo)~ (VO) - J (vo)~ (VO) P P P P 

(A-17) 

(A-18) 

2k 

6 vo = - 

and the primes denote the derivative with respect to v. For the 

solutions in this paper, P=k with kD=l. Once the B's are 

determined, the electric field at any point in the layer is 

computed by using equation (16). If p is not zero or a positive 

integer, the general solution of equation (A-14) becomes 

Ex(v) = D J (v) + D J- (v) (A-19) 1P 2P 

where the constants D1 and D2 are determined in the manner just 

described. 

15 
REFERENCES 

1. Stratton,J.A.; Electromagnetic Theory, McGraw-Hill Book 

Company, New York, N.Y., 1941. 

2. Barrar,R.B.; and Redheffer,R.M.; On Nonuniform Dielectric 

Media," IRE Transaction on Antennas and Propagation, 

Vol.AP-3,July 1955. 

3. Budden,K.G.; Radio Waves in the Ionosphere, Cambridge 

University Press, 1961. 

4. Wait,J.R.; Electromagnetic Waves in Stratified Media, Pergamon 

Press, New York,N.Y.,1970. 

5. Sharpe,C.B.; The Scattering Approach to the Synthesis of 

Nonuniform Lines,l# Research Institute, University of 

Michigan, Ann Arbor, Rept. No.119, March 1961. 

6. Abramowitz,M.: and Stegun,I.A.;editors; Handbook of 

Mathematical Functions, National Bureau of Standards, U.S. 

Government Printing Office, Washington,D.C.,May 1968. 

7. Richmond,J.H.; Transmission Through Inhomogeneous Plane 

Layers," IRE Transactions on Antennas and Propagation, May 

1962. 

16 
Y 

t 

Region 

ZI Region Ill 

. Figure 1 .-Perpendicularly-polarized plane wave incident 

on an inhomogeneous plasma layer. 

17 
I 

I 

Figure 2.-Linear 

I 

I 

I z< 

D 

Figure 3.-Exponen 

2 

-1 

prof i I e. 

\ 

e 

1 

:ial profile. 

. 

18 
Report Documentation Page 

1. Report No. 

NASA TM-101609 

2. Government Accession No. 

9. Performing Organization Name and Address 

I 

NASA Langley Research Center 

Hampton, VA 23665-5225 

unclassified 

12. Sponsoring Agency Name and Address 

19 A0 3 uncl ass if ied 

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

15. Supplementary Notes 

i 

3. Recipient's Catalog No. 

i 5. Report Date 

August 1989 

6. Performing Organization Code 

8. Performing Organization Report No. 

10. Work Unit No. 

583-01-11-12 

11. Contract or Grant No. 

13. Type of Report and Period Covered 

Techni cal Memorandum 

14. Sponsoring Agency Code 

16. Abstract 

Numerical solutions of the differential equation which describe the electric field 

within an inhomogeneous layer of permittivity, upon which a perpendicularlypolarized plane wave is incident, are considered. Richmond's method and the 

Runge-Kutta method are compared for linear and exponential profiles of permittivities. These two approximate solutions are also compared with the exact solutions. 

17. Key Words (Suggested by Authork)) 

di fferenti a1 equation 

n urne r i ca 1 

Runge-Kutta 

18.Distribution Statement 

unclassified - unlimited 

I Subject Category 32