Comp081011.qxd


The Journal of Engineering Research  Vol. 7, No. 1 (2010)  62-69

1.  Introduction

The most popular techniques used for numerical sim-
ulations of various electronic devices are the Finite
Difference and the Finite Element Methods (FDM and
FEM). Although, some other techniques, such as the
Boundary Integral and Monte Carlo methods have been
used in some cases, the two former methods are the most
widely used despite of some of their disadvantages. 

The main drawbacks of these methods are long run-
ning  time, large memory space requirements and conver-
____________________________________
Corresponding author’s e-mail: hadj@squ.edu.om

gence problems. In the present work, an alternative semi-
numerical approach is used for the simulation of electron-
ic planar devices. The approach is based on a mathemati-
cal method known as the Method of Lines (MoL) (Ames,
1977, Sadiku, et al.  2000; Vietzorreck, et al.  2000;
Barcz, et al.  2003; Vietzorreck, et al. 2004; Pascher, et al.
2005; Yan, et al.  2005; Gonzalo, et al.  2006; Plaza, et al.
2006 and  Chen, et al.  2007)  which is a differential-dif-
ference technique. It is a versatile technique developed by

Determination of Potential Profile In Planar Electronic
Structures Using A Semi-Analytical Technique

Hadj Bourdoucen**a,  Mokrane Dehmasb and El-Bachir Yallaouic

*a Department of Electrical and Computer Engineering, College of Engineering, Sultan Qaboos University, PO Box 33, PC 123,
Al-Khoud, Sultanate of Oman

**a Communication and Information Research Centre, Sultan Qaboos University, Muscat, Sultanate of Oman
b University of Boumerdes, Department of Electrical Engineering,  35000, Algeria

c Department of  Mathematics and Statistics, College of Science, Sultan Qaboos University, PO Box 36, PC 123, 
Al-Khoud, Sultanate of Oman

Received  11 October 2008; accepted  13 September 2009

Abstract:  In this paper, a semi-analytical technique known as the Method of Lines (MoL) with uniform and
non-uniform discretization schemes is developed. The aim is to determine static potential profile in planar elec-
tronic structures. Even though this method has been known for some time, there has been reported work on its
application to planar semiconductor device analysis for voltage profile determination.  Since most current elec-
tronic devices are manufactured using planar and quasi-planar technology, the proposed algorithm is well suit-
ed for device analysis prior to fabrication. Compared with known popular methods such as Finite Difference and
Finite Elements methods, the proposed technique is relatively simple, more accurate and unlike other methods,
has no convergence problem. In addition to this, its semi-analytical nature, which consists of reducing one com-
puting dimension, allows saving significant memory and computation time. Typical planar electronic structures
are considered to demonstrate their suitability for these devices, and the obtained results are presented and dis-
cussed.

Keywords:  Planar structure, Method of lines, Potential profile, MoL

 á«∏«∏ëàdG ¬Ñ°T á«æ≤àdG ΩG~îà°SEÉH ájƒà°ùe á«fhÎμdEG πcÉ«¡d ~¡÷G íeÓe ~j~–

 …hÓ©j Ò°ûÑdG , ¢SɪMO ¿Gô≤e ,ø°ShOQƒH êÉM

::áá°°UUÓÓÿÿGG) •ƒ£ÿG á≤jôW º°SEÉH áahô©ŸG h á«∏«∏ëàdG ¬Ñ°T á«æ≤àdG ôjƒ£J ” åëÑdG Gòg ‘MoL~¡÷G íeÓe ~j~– ¢Vô¨H ∂dPh ᪶àæe ÒZh ᪶àæe á∏°üØæe äÉ££fl ™e (

~©jh .ájƒà°ùe á∏°Uƒe áÑ°T §FÉÑæd ~¡÷G íeÓe ~j~– h π«∏ëàd IQƒ°ûæe çÉëHCG ~Lƒj ’ áfCÉa ≥HÉ°S âbh òæe áahô©e á≤jô£dG √òg ¿CG ºZôHh . ájƒà°ùe á«fhÎμdEG πcÉ«¡d øcÉ°ùdG

¥ô£dÉH á≤jô£dG √òg áfQÉ≤à h .ájƒà°ùe áÑ°T hCG ájƒà°ùe äÉ«æ≤J ΩG~îà°SÉH ™æ°üJ §FÉÑædG ∂∏J Ö∏ZG ¿G å«M É¡©«æ°üJ πÑb á«fhÎμdC’G §FÉÑædG á°SGQ~d ÉÑ°SÉæe ìÎ≤ŸG ΩRQGƒÿG

.iôNC’G ¥ô£dG πãe ÜQÉ≤àdG á∏μ°ûe øe ≈fÉ©J ’ h ábO ÌcCG É¡fG ɪc É«Ñ°ùf §°ùHCG áMÎ≤ŸG á«æ≤àdG ¿CG ~‚ IOh~ÙG ô°UÉæ©dG á≤jôW h IOh~ÙG ¥QGƒØdG á≤jôW πãe á©FÉ°ûdG

áMÎ≤ŸG á≤jô£dG ≥«Ñ£J ”h .ÜÉ°ù◊G øeRh IôcGP øe πμd ≠dÉH ¢ ØîH íª°ùj ɇ ~MGh ≈Ñ°SƒM ~©H ¢ «ØîJ πª°ûJ á«∏«∏– áÑ°T á©«ÑW É¡d áMÎ≤ŸG á«æ≤àdG ¿EÉa ≥Ñ°SÉŸ áaÉ°VC’ÉH

 .É¡«∏Y ∫ƒ°ü◊G ” ≈àdG íFÉàædG á°ûbÉæeh í«°VƒJ ”h »°SÉ«b iƒà°ùe ≈fhÎμdCG πμ«g

áá««MMÉÉààØØŸŸGG  ääGGOOôôØØŸŸGG.~¡÷G íeÓe ,•ƒ£ÿG á≤jôW ,ájƒà°ùŸG πcÉ«¡dG :



63

The Journal of Engineering Research  Vol. 7, No. 1 (2010)  62-72

R. Pregla and co-workers for the analysis of planar
microwave structures (Schultz, 1980, Pregla and Pascher,
1989, Pregla, 2002).  For a given partial differential equa-
tion subject to some boundary conditions, all but one of
the independent variables are discretized. This allows
finding a one-dimensional analytical solution. The semi-
analytical nature of the MoL makes the computational
efforts much less intensive than the above-cited methods
applied so far to microwave structures. This method has
been used recently for the necessity of reducing long cal-
culation times and therefore to  significantly reduce the
computation time (Zhou, et al.  2000; Pascher, et al.  2005;
Yan, et al.  2005; Gonzalo, et al.  2006; Plaza, et al.  2006
and Chen, et al.  2007). The method of lines has been used
by several authors for the analysis of microwave struc-
tures in both quasi-static and full-wave modes (Schultz,
1980; Pregla and Pascher, 1989; Pregla, 2002; Barcz, et
al.  2003; Vietzorreck and Pascher, 2004 and Yan, et al.
2005).  However, the method has not been extended to the
voltage profile determination of planar semiconductor
devices. In the following sections, the application of the
MoL to the analysis of planar structures will be demon-
strated based on a two-dimensional Poisson's equation.
Since the discretization schemes are important parameters
in the efficiency of the algorithm, the uniform and non-
uniform schemes will both be presented.  Finally, some of
the obtained results for basic semiconductor planar struc-
tures will be depicted.

2.  Technique of Analysis

Let us introduce the approach of the Method of Lines
based on solving the very popular Poisson's equation
given by Eq.  (1) for the domain shown in Fig. 1.

(1)

If we uniformly discretize the variable x, then the func-
tions (x ,y) and  f(x ,y) in Eq. (1) are transformed into the

sets (xi,y) and  f(xi,y) along the lines  x=i.h where
i=1,2,3…N. 

N is the total number of lines within the structure and
h is the discretizing size interval as shown in Fig. 1.

This operation transforms Eq. (1) into a system of N
differential equations of the form:

(2)

If we let (xi,y) = i(y) and  f(xi ,y) = fi(y) then, the ith
difference approximation on the discretized variable may
be written as:

(3.a)

and its second derivative may be written as:

(3.b)

If we express the vectors of functions  and F as: 

(4.a)

and

(4.b)

then, Eq. (2) can be written in the following matrix
form:

(5)

where P is a second order differential operator, which
is an NxN tri-diagonal matrix of the form:

x i  
x i1

  


2

2
1

x
i 

 
 


2

2x
i

 

 
x i

 
x

i 1

 


2

2
1

x
i

i1

h  h/2. . .

y

i1

. . .

Figure 1.  Equidistant discretization pattern using the method of lines

     
2 2

2 2
x, y x, y

f x, y
x y

   

 
 

     
2 2

2 2
i i

i
x , y x , y

f x , y
x y

   

 
 

1i i

ix h
 


 

2
1 1 1

2 2
2i i i i i

i

x x
hx h

 
    


  


 

 

 1 2
t

N, ,...,   

 1 2
t

NF f , f ,..., f

2

2 2
1d

P F
dy h


 



64

The Journal of Engineering Research  Vol. 7, No. 1 (2010)  62-69

(6)

where p1 and p2 are determined by the left and right
boundaries of the structure of Fig. 1.

Since P is a real-symmetric and tri-diagonal matrix,
The following two possible conditions on  0 and N+1
may be met:

there exists a nonsingular NxN matrix T such that the
matrix   given by:

(7)
is diagonal and where the elements  i are  the eigenval-
ues of  P. 

T is a matrix of the corresponding eigenvectors and Tt
is its transpose. 

For a non-uniform discretization scheme, the approach
presented above has to be slightly modified. 

Hence, the second derivative expressed by Eq. (3.b) for
the uniform discretization will now be expressed as given
below for the non-uniform case:

(8)

where ei = (h1 + hi-t) / 2 denotes the ith interval size
between the dotted lines shown in Fig. 2 and  hi is the
non-uniform discretization interval.

(9)

where h is an assumed normalized discretization interval.
Considering the Dirichlet-Dirichlet lateral boundary

conditions of the analyzed structure, which means  
 , Eq. (9) takes the following matrix form

(10)

where

(11)

(12)

and

(13)

whereas the matrix D of dimension (N+1xN) is a first
order difference operator given by:

(14)

(15)
which takes the matrix form:

(16)

where

(17)

and

(18)

1

2

1 0
1 2 1

1 2 1
0 1

p

P ... ... ...

p

 
 
 
 
  

1. either 0  and/or 01N   where in this case 

1p  and/or 2p2   (Dirichlet condition): 

2. or 
0x


 and/or 0

x N





 where 1p  and/or 

1p2    (Neumann condition). 

tT PT 

2
1

2
1 2i i

ii

x x
 ,      i , , .... N .

ex

 
  





 

ei ei+1

hi-1 hi 

 i i 1


x i1


x i1


x i

i1

x

y

Figure 2.  Non-equidistant discretization pattern with
the method of lines

     The difference Eq. (3.a) for the non-uniform case is 
firstly normalized by multiplying both of its sides by 

hhh i  and to get the relation: 

 1i i i
ii

h h
h

h x h


 
 

 
    

 

1
h x hhr r D 
 

0 1h
i

h
r diag      ;      i , ,... N .

h
 

   
 

 1 2
t

N, , k ,    

0 1

t

x
N

, , ..., 
x x x

  


  
 

   
 

1 0
1 1

1
1

0 1

D ...
...

 
 
  
 
  

      Similarly, the second order operator  given by Eq. 
(8) is  normalized by multiplying both of its sides by 
the factor ih e / h  to get the following expression:  

2

2
1

1 2i
i i ii

e h
h     ;    i , , ..., N .

h e x xx
   

  

   
          

1 t
e xx e xhr r D 
  

1e
i

h
r diag   ,        i , 2, ...N .

e
 

   
 

2 2 2

2 2 2
1 2

t

xx
N

, , ..., 
x x x

     


  

 
 
 
 



65

The Journal of Engineering Research  Vol. 7, No. 1 (2010)  62-69

and Dt is the transpose of D.
Inserting the expression of x given by Eq. (10) into

Eq. (16) gives:

(19)

If we let

(20)

and
(21)

then Eq. (19) may be written as:

(22)

Inserting the expression of  xx from (19) into the sys-
tem given by Eq. (2) gives the following matrix form

(23)

where f is a vector of functions obtained  from discretiz-
ing  f (x,y).

(24)

where
(25)

and

(26)

Since Dxx is a tri-diagonal matrix, then the system
given by Eq. (24) needs to be decoupled before being
solved.  Furthermore, since Dxx is real-symmetric,  there
exists an orthogonal matrix T such that the elements of the
diagonal matrix   are the eigenvalues of P.

T in this case, is also the matrix of the corresponding
eigenvectors and Tt its transpose. Using the above result,
the system given by Eq. (24) is decoupled by pre-multi-
plying it by Tt and obtaining the following system:

(27)

or simply

(28)

where

(29)

and

(30)

(31)

The general solution of each equation of the system
(31) takes the form:

(32)

(33)

where the coefficients ai are determined by substituting
Vi of  Eq. (31) by the above expression.

In the end, the vector of potentials  is obtained from
V using the following expression:

(34)

which expresses analytically the solutions along the
discretization lines as a function of  y.

Based on the above algorithm, a computer program has
been developed for the simulation of semiconductor
devices to determine the potential distribution within pla-
nar structures. It includes the following steps:

1.  Enter the structure geometry and the dimension of each
layer, the type of boundary conditions, the electric
charge profile as well as the discretization interval h.

2.  Perform discretization by finding the total number of
lines and construct the second order operator P.

3.  Determine the eigenvalues of the matrix P, the ele-
ments of   and the corresponding matrix of eigenvec-
tors T.

4.  Find in the transformed domain the space charge den-
sity and homogeneous horizontal surface boundary
vectors by pre-multiplying by  Tt.

5.  Compute the constants Ai , Bi by considering the field
and potential interface continuity conditions as well as
the surface boundary conditions.

6.  Perform the inverse transformation to find the analytic
solution in the original domain along each discretiza-
tion line.

   2 1 1te xx e h h e eh r r D r r Dr r   

x h eD r Dr

t
xx x xD D D 

 2 1 1e xx xx eh r D r  

 
2

1
2 2

1
e xx e

d
r D r f

dy h


 

      Multiplying both sides of Eq.  (23) by re
1 yields:  

2

2 2
1

xx
d

D G
dy h


 

1
er 


1
eG r f


    
2

2 2
1

t
t t t

d T
T PT T T G

dy h


 

2

2 2
1d V

V F
dy h

 

tF T G

tV T 

     Assuming  i i 
2  then Eq.  (28) is further 

transformed into a set of N ordinary differential 
equations as 

22

2
1i i i i

d V
V F i ,..., N

hdy

 
   
 

i i
i i i PiV A cosh y B sinh y Vh h

    
     

   

where Ai and Bi are constants depending on the 
horizontal side boundaries. The specific solution VPi  is 
expressed as a linear combination of the functions Fi  
and its linearly independent derivatives 

 n
iF  as 

     1 2
0 1 2

m
Pi i i i i i i i imV a F a F a F ... a F ...     

e er  r  T  V  



66

The Journal of Engineering Research  Vol. 7, No. 1 (2010)  62-69

3.  Case Studies

In this section, the developed algorithm is applied to
a set of typical planar structures used in current electron-
ics technology. The objective of the simulation is to deter-
mine the potential profile throughout the structures having
different physical, geometrical and boundary conditions.

Using a uniform discretization scheme with an interval
size h = 1  m  results in a total of 30 intervals,  the fol-
lowing lateral boundary conditions have been considered: 

a) Neumann-Dirichlet
b) Dirichlet-Neumann
c) Dirichlet-Dirichlet
d) Neumann-Neumann 

The results showing the obtained potential profile con-
tours are depicted in Fig. 4.

With reference to this figure, one can observe that the
change in the potential profile is strongly related to the
boundary conditions. This is in accordance with most
practical situations where a boundary of a semiconductor

 

H2

H1

l2
L

Va (applied voltage) 
 y

x

(N)

Layer 1
1 

Layer 2
2 

l1

(N)

Figure 3.  Cross sectional view of a two-layer planar semiconductor structure.  (Dimensions are given in text)

0 5 10 15 20 25 30
0

5

10

15

 x (µm)

  y (µm)

 a) N-D

19 v

16 v

13 v

10 v

7 v

4 v

1 v

 
 

 0 5 10 15 20 25 30
 0

 5

 10

15

 x (µm)

 y (µm)

 b) D-N

 19 v

16 v

13 v

10 v

7 v

4 v

1 v

 

0 5 10 15 20 25 30
0

5

10

15

 x (µm)

  y (µm)

 c) D-D

19 v

16 v

13 v

10 v

7 v

 4 v

1 v

 

 0 5 10 15 20 25 30
 0

5

10

15

 x (µm)

 y (µm)

d) N-N

19 v

16 v

13 v

10 v

7 v

4 v

1 v

 
 

 Figure 4.  Potential profiles in the structure in Figure 3 with Neumann-Dirichlet,  b) Dirichlet-Neumann, 
c) Dirichlet-Dirichlet and  d) Neumann-Neumann lateral boundary conditions

      The effect of boundary conditions on the potential 
profile has been first considered by finding the potential 
profile throughout a two-layer planar structure as 
shown in Fig. 3.  The electrical, physical and 
geometrical constants of this structure are: 1 = 10

-11 
C/cm3;  2 = -2 x 10

-11 C/cm3; Va = 20 V; 1 =  2 = 119 
0; L = 30 m; 11 = 8 m; 12 = 10 m; H1 = 10 m; H2 = 
5 m and the value of 0 is 8.85x10-12 F/m. 



67

The Journal of Engineering Research  Vol. 7, No. 1 (2010)  62-69

device is either forced to have a specified voltage or left
free.

After this, two-layer symmetrical structures shown in
Fig. 5 have been simulated. Note that in both structures
the layers near to the groundsides have a constant charge
density , whereas the charge profile throughout the second
layers varies linearly with y as   A(y-H1) + B.  The
dimensions and the physical parameters for both struc-
tures are:  1 = 3.5 x 1011 C / cm3; A = -1013 C / cm4 B =
-1011 C / cm3; = 11.0  ; L = 30 m mH2 =
15 m; and Va = 20 V.   For the structure 5-a,  11 = 12 =
10 m and the boundary conditions are Dirichlet-
Dirichlet; whereas for the structure 5-b, 11 = 12 = 6 m
and the boundary conditions are Neumann-Neumann.

The analysis is carried out with 63 uniform intervals

for the device 5-a and 65 intervals for the device 5-b.
Figure 6 shows the simulation results of the obtained
potential profiles for both structures 5-a and 5-b. 

Due to symmetry, only half of these structures need to
be simulated by considering the Neumann condition at the
axis of symmetry. It has been verified that the analytical
results are the same for all symmetric lines with respect to
the axis of symmetry. Hence, for symmetrical devices the
numerical effort can be significantly reduced by reducing
the simulated space and therefore, reducing the total num-
ber of lines by a factor of two.

A two-dimensional metal oxide semiconductor (MOS)
capacitor has been also simulated using the developed
algorithm. The structure consists of a metallic gate, a P-
type silicon substrate and an oxide film layer as shown in
Fig. 7.

 

 

 

 

x (a) 

H2 

H1 

l1 
L 

(N) 

1  

2  l2 

x 

y 

H2 

H1 

l1 
L 

Va 

1  

2  l2 

 (b) 
Figure 5.  Cross sectional view of two-layer semiconductor symmetrical structures

0 5 10 15 20 25 30 
0  

5.0 

10 

15 

x (µm ) 

y (µm) 

 (a) 

16 v 

13 v 

9 v 

5 v 

1 v 

 

 

0 5 10 15 20 25 30 
0 

5 

10 

15 

x (µm) 

 y (µm) 

 (b) 

16 v 16 v 

12 

8 v 

4 v 

1 v 

 
 

(a)                                                                                 (b)

Figure 6.  Potential profiles for the structures of Figures 5-a  and  5-b  respectively



68

The Journal of Engineering Research  Vol. 7, No. 1 (2010)  62-69

Because of the presence of a layer of silicon dioxide
layer (which is a dielectric material) between the gate and
the substrate, this device exhibits the properties of a
capacitor. It is assumed in this analysis that the metal and
the semiconductor work functions are equal, and that there
is no electric field either within the oxide or at the oxide-
silicon interface. Hence, the electric field is supposed to
be zero whenever the applied voltage is zero.

Since the two-dimensional model of the considered
device is symmetric, the investigation is limited to the left
part of the device as illustrated in Fig. 7.  When a positive
voltage is applied at the electrode, free majority carriers
(holes) are repelled by the induced electric field.
Equilibrium condition requires that the potential is zero at
the self-adjusting boundaries of the resulting fully-deplet-
ed layer and elsewhere in the neutral region. The lateral
boundary conditions of the silicon and oxide layers are of
Neumann type.

The discretization pattern is formed by 10 concentrat-

ed equidistant lines crossing the curved area and 31 large-
ly spaced lines elsewhere. The space charge layer is sub-
divided into 10 sub-layers to perform a stair-step approx-
imation on the curved boundary. The dimensions used are:
L = 35 m; l = 17 m; 0.1 < Xo < 5.0 m; 30 < Xs < 50
m; 1014 < NA < 1018 cm-3; 10 < VG < 100 (VG <
Vbreakdown ).

The program has been executed for different  Xo and
NA and the results obtained are shown in Fig. 8. The
results shown in this figure agree with both the previous-
ly stated boundaries and, the physical requirements. An
important observation is that considering only a part of the
neutral region can still save a lot of numerical efforts. This
is because the potential within the whole neutral region is
zero. This last point is particularly important when the left
lateral side is far from the space charge region boundary,
and consequently, the number of discretization lines can
also be reduced.

 VG

x

y

Xo

Xs

Xd

P-

SiO2

Depletion Neutral 

l
L

Figure 7.  Two-dimensional model of the MOS capacitor used in the present analysis

0 
10 

20 
30 

40 

0 

5 

10 

15 
-50 

0 

50 

100 

150 

x (μm) y (μm) 

V 

Figure 8.  Three-dimensional representation of the potential distribution in the device of Figure 7



69

The Journal of Engineering Research  Vol. 7, No. 1 (2010)  62-69

4.  Conclusions

The Method of Lines (MoL) is an efficient numerical
tool for solving partial differential equations. Its strength
lies in its assured convergence and its semi-analytical
nature for both uniform and non-uniform discretization
schemes. As it appears from the analysis done in this
work, the MoL can be applied to electronic devices to find
static potential profiles for different structures. A direct
application of the MoL is to determine the breakdown
voltages of power electronic devices needed for design
prior to manufacturing.  The developed algorithm allows
significant saving in memory storage and computation
time with assured convergence. This facilitates its exten-
sion to more sophisticated electronic structures  as well as
its integration into current computer-aided design tools. 

References

Ames, W.F., 1977, "Numerical Methods for Partial
Differential Equations," Academic Press, New York.

Barcz, A., Helfert, S.F. and  Pregla, R., 2003, "The
Method of Lines Applied to Numerical Simulation of
2D and 3D Bandgap Structures," Proceedings of  5th
International Conference on Transparent Optical
Networks, 1,  pp. 126-129. 

Chen, R. S.,  Fang, D.G. and Li, X.G., 2007,  "Analysis of
Open Microstrip Lines by MOL," International
Journal of Microwave and Millimeter-Wave
Computer-Aided Engineering, Vol. 3(2), pp. 109 -
113.

Gonzalo, P., Ricardo, M. and Francisco, M., 2006, "Quasi-
TM MoL/MoM Approach for Computing the
Transmission-Line Parameters of Lossy Lines," IEEE
Trans. Microw. Theory Technology,  Vol. 54(1), pp.
198-209.

Pascher, W., Pregla, R. and Vietzorreck, L., 2005, "Fast
Full-Wave Analysis of Distributed MEMS
Transmission Lines by the MoL," International
Conference on Wireless Communications and

Applied Computational Electromagnetics,
IEEE/ACES, 3-7 ,  pp. 763 - 766.

Plaza, G., Marques, R. and Medina, F., 2006, "Quasi-TM
MoL/MoM Approach for Computing the
Transmission-Line Parameters of Lossy Lines," IEEE
Transmission Microw. Theory Tech.nology, Vol.
54(1), pp. 198-209.

Pregla, R.  and Pascher, W., 1989, " Numerical Technique
for Microwave and Millimeter Passive Structures,"
Ed. John Willey and Sons. 

Pregla, R., 2002, "Efficient and Accurate Modeling of
Planar Anisotropic Microwave Structures by the
Method of Lines,"  IEEE Transmission Microwave
Theory and Techniques, Vol. 50, pp.  1469-1479.

Sadiku, M.N.O. and Garcia, R.C., 2000, "Method of Lines
Solution of Axisymmetric Problems," Southeastcon
2000. Proceedings of the IEEE, pp. 527 - 530.

Schultz, U., 1980, "The Method Of Lines - A New
Technique for the Analysis of Planar Microwave
Structures," PhD Thesis, Fern Univ. Hagen, Germany. 

Vietzorreck, L., Coccetti, F., Chtchekatourov, V.  and
Russer, P., 2000, "Numerical Methods for the High-
Frequency Analysis of MEMS Capacitive Switches,"
Proceedings of the II. Topical Meeting on Silicon
Monolithic Integrated Circuits in RF Systems 26 - 28.
April 2000 Garmisch Germany, pp. 123-124. 

Vietzorreck, L. and  Pascher, W., 2004, "Efficiency
Enhancement by Reduction of Modal Complexity in
the Analysis of Cascaded Planar Circuits by the
MoL," Microwave Symposium Digest, IEEE MTT-S
International, 3,  pp- 1423-1426.

Yan, L.,  Hong, W., Wu, K. and Cui, T.J., 2005,
"Investigations on the Propagation Characteristics of
the Substrate Integrated Waveguide Based on the
Method of Lines," Microwaves, Antennas and
Propagation, IEE Proceedings, Vol. 152(1),   pp. 35 -
42. 

Zhou, G.R., Fang, D.G. and  Feng, N.N., 2000, "Entire
Domain Based Diakoptic Method of Lines," CEEM
2000 Proceedings, Asia-Pacific Conference on
Environmental Electromagnetics, pp. 270-273.