Review: Numerical Methods for Differential Equations

ECE Generic and 2331

Han Q. Le(c)

Segment 2 - System and higher-order equations

1. Basic idea: multifunction and variable

More often than not, the real world problems involve many parameters. Instead of just a function y[x], we have u[x], v[x], etc... that are related to each other.

How do we apply techniques such as Euler's or R-K's to these problems? It is actually simpler than we think. For convenience, we treat the functions as a single entity, a vector function, like this:
((d (Subscript[y, 1]))/(d x))=(Subscript[f, 1])(x, (Subscript[y, 1]), (Subscript[y, 2]),... ,(Subscript[y, n]))
((d (Subscript[y, 2]))/(d x))=(Subscript[f, 2])(x, (Subscript[y, 1]), (Subscript[y, 2]),... ,(Subscript[y, n]))
...
((d (Subscript[y, n]))/dx)=(Subscript[f, n])(x, (Subscript[y, 1]), (Subscript[y, 2]),... ,(Subscript[y, n]))

or:
(d/(dx))(({{(Subscript[y, 1])}, {(Subscript[y, 2])}, {.}, {(Subscript[y, n])}}))=(({{(Subscript[f, 1])(x, (Subscript[y, 1]), (Subscript[y, 2]),... ,(Subscript[y, n]))}, {(Subscript[f, 2])(x, (Subscript[y, 1]), (Subscript[y, 2]),... ,(Subscript[y, n]))}, {.}, {(Subscript[f, n])(x, (Subscript[y, 1]), (Subscript[y, 2]),... ,(Subscript[y, n]))}}))  
or:
(d/(dx))Y=F(x,Y) 
where Y is a vector and so is F(x,Y). 
In some field, this is also known as a state-space model.
So far, we haven’t done anything new, just a convenient notation.
But we can apply all previous methods to this problem

2. Euler's method for multifunction

2.1 Basic idea

Here is our equation: (d/(dx))Y=F(x,Y) 
For illustration, let's take Y to be of 3-dimension.
From Taylor's series expansion, we recall that:
(({{(Subscript[y, 1])((Subscript[x, i+1]))}, {(Subscript[y, 2])((Subscript[x, i+1]))}, {(Subscript[y, 3])((Subscript[x, i+1]))}}))  ≈ (({{(Subscript[y, 1])((Subscript[x, i]))}, {(Subscript[y, 2])((Subscript[x, i]))}, {(Subscript[y, 3])((Subscript[x, i]))}})) + ((Subscript[x, i+1]) - (Subscript[x, i])) (d/(dx))(({{(Subscript[y, 1])}, {(Subscript[y, 2])}, {(Subscript[y, 3])}}))
So, we estimate:
(({{(Subscript[y, 1])((Subscript[x, i+1]))}, {(Subscript[y, 2])((Subscript[x, i+1]))}, {(Subscript[y, 3])((Subscript[x, i+1]))}}))  ≈ (({{(Subscript[y, 1])((Subscript[x, i]))}, {(Subscript[y, 2])((Subscript[x, i]))}, {(Subscript[y, 3])((Subscript[x, i]))}})) + h (({{(Subscript[f, 1])((Subscript[x, i]);(Subscript[y, 1])((Subscript[x, i])),(Subscript[y, 2])((Subscript[x, i])),(Subscript[y, 3])((Subscript[x, i]))) }, {(Subscript[f, 2])((Subscript[x, i]);(Subscript[y, 1])((Subscript[x, i])),(Subscript[y, 2])((Subscript[x, i])),(Subscript[y, 3])((Subscript[x, i])))}, {(Subscript[f, 3])((Subscript[x, i]);(Subscript[y, 1])((Subscript[x, i])),(Subscript[y, 2])((Subscript[x, i])),(Subscript[y, 3])((Subscript[x, i])))}}))
We see that apart for writing everything 3 times, everything else is the same

2.2 Example

(d/(dx))(({{(Subscript[y, 1])}, {(Subscript[y, 2])}, {(Subscript[y, 3])}})) = (({{(Subscript[f, 1])(x, (Subscript[y, 1]), (Subscript[y, 2]),... ,(Subscript[y, n]))}, {(Subscript[f, 2])(x, (Subscript[y, 1]), (Subscript[y, 2]),... ,(Subscript[y, n]))}, {(Subscript[f, 3])(x, (Subscript[y, 1]), (Subscript[y, 2]),... ,(Subscript[y, n]))}}))= (({{-3((Subscript[y, 2])-(Subscript[y, 3]))}, {15(Subscript[y, 1])-(Subscript[y, 1]) (Subscript[y, 3])}, {1.5(Subscript[y, 1]) (Subscript[y, 2])-(Subscript[y, 3])}}))
If at x=0; (({{(Subscript[y, 1,1])}, {(Subscript[y, 2,1])}, {(Subscript[y, 3,1])}})) =(({{1}, {-1}, {1}})) ; then at x=h = 0.1, we have an estimate:
(({{(Subscript[y, 1 new])}, {(Subscript[y, 2 new])}, {(Subscript[y, 3 new])}})) ==  (({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}, {(Subscript[y, 3 old])}})) + h * (({{-3((Subscript[y, 2 old])-(Subscript[y, 3 old]))}, {15(Subscript[y, 1 old])-(Subscript[y, 1 old]) (Subscript[y, 3 old])}, {1.5(Subscript[y, 1 old]) (Subscript[y, 2 old])-(Subscript[y, 3 old])}}))

2.2.1 Algorithm

Manipulate[
If[!TrueQ[{xendold,hdBold}=={xend,hdB}],
h=Max[10^(0.1 hdB),0.001 xend] ;
If[!TrueQ[hold==h],
np=Ceiling[xend/h];
If[np≤ 4, h=xend/4.;np=Ceiling[xend/h];,];
(({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}, {(Subscript[y, 3 old])}}))  = (({{1}, {1}, {0}})); x = 0. ; 
{yvec,xarr}=Reap[For[i=1, i≤ np, {
(({{(Subscript[y, 1 new])}, {(Subscript[y, 2 new])}, {(Subscript[y, 3 new])}})) =  (({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}, {(Subscript[y, 3 old])}})) + h *(({{-3((Subscript[y, 2 old])-(Subscript[y, 3 old]))}, {15(Subscript[y, 1 old])-(Subscript[y, 1 old]) (Subscript[y, 3 old])}, {1.5(Subscript[y, 1 old]) (Subscript[y, 2 old])-(Subscript[y, 3 old])}})) ;
x  = x +h ;
Sow[{(Subscript[y, 1 new]),(Subscript[y, 2 new]),(Subscript[y, 3 new])},A];Sow[x,B];
		   (({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}, {(Subscript[y, 3 old])}}))  = (({{(Subscript[y, 1 new])}, {(Subscript[y, 2 new])}, {(Subscript[y, 3 new])}}))  },i++]][[2]];
PrependTo[yvec,{1,1,0}];
yvecgr=Table[{Hue[((i-0)/(Length[yvec]))1,0.8],Line[{yvec[[i]],yvec[[i+1]]}]},{i,Length[yvec]-1}];
hold=h,];
{xendold,hdBold}={xend,hdB},];

Graphics3D[{Thick,yvecgr}]

,Grid[{
{Control[{{hdB,-20,Button["step size (dB)",BaseStyle-> {14, Bold,Black,FontFamily-> "Arial"}
,Background-> Hue[0.16,0.7,1]]
},-40,-20,ImageSize-> Medium}]
,Control[{{xend,10,Button["x end point",BaseStyle-> {14, Bold,Black,FontFamily-> "Arial"}
,Background-> Hue[0.04,0.7,1]]},5,10.,ImageSize-> Medium}]
}}]
,Initialization:>{
{xendold,hdBold}={-1.,0}

}]

3. Example of Runge-Kutta

3.1 Problem set up

Now we see that besides the fact that we have to keep track of more than one functions, the rest is essentially the same as what we did before. Here we'll look at a 4th order R-K for a 2-function case:
(d/(d t))Lynx = -0.9*Lynx + 0.5*Lynx*Hare
(d/(d t))Hare = 1.1*Hare - 0.7*Lynx*Hare

3.2 Solution

We remember that the fourth order expression in R-K method is:
(Subscript[y, i+1]) = (Subscript[y, i]) + (1/6) ((Subscript[k, 1]) + 2 (Subscript[k, 2]) +2 (Subscript[k, 3]) +(Subscript[k, 4]) ) h ;
where:
(Subscript[k, 1]) =  f((Subscript[x, i]),(Subscript[y, i])) ;
(Subscript[k, 2]) =  f((Subscript[x, i]) + (1/2)h,(Subscript[y, i]) +  (1/2)(Subscript[k, 1]) h) ;
(Subscript[k, 3]) =  f((Subscript[x, i]) +(1/2) h,(Subscript[y, i]) + (1/2)(Subscript[k, 2]) h  ) ;
(Subscript[k, 4]) =  f((Subscript[x, i]) + h,(Subscript[y, i]) + (Subscript[k, 3]) h  ) ;
We set up now for a vector
(({{(Subscript[y, 1 new])}, {(Subscript[y, 2 new])}}))=(({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}}))+(1/6) h (({{(Subscript[k, 1,1])+2 (Subscript[k, 1,2])+2 (Subscript[k, 1,3])+(Subscript[k, 1,4])}, {(Subscript[k, 2,1])+2 (Subscript[k, 2,2])+2 (Subscript[k, 2,3])+(Subscript[k, 2,4])}}));
where:
(({{(Subscript[k, 1,1])}, {(Subscript[k, 1,2])}, {(Subscript[k, 1,3])}, {(Subscript[k, 1,4])}})) = (({{(Subscript[f, 1])((Subscript[x, i]),(Subscript[y, 1,i]), (Subscript[y, 2,i]))}, {(Subscript[f, 1])((Subscript[x, i]) + (1/2)h,(Subscript[y, 1,i]) +  (1/2)(Subscript[k, 1,1]) h , (Subscript[y, 2,i]) +  (1/2)(Subscript[k, 2,1]) h)}, {(Subscript[f, 1])((Subscript[x, i]) +(1/2) h,(Subscript[y, 1,i]) + (1/2)(Subscript[k, 1,2]) h  , (Subscript[y, 2,i]) +  (1/2)(Subscript[k, 2,2]) h)}, {(Subscript[f, 1])((Subscript[x, i]) + h,(Subscript[y, 1,i]) + (Subscript[k, 1,3]) h  , (Subscript[y, 2,i]) + (Subscript[k, 2,3]) h)}})) ; and similarly:
(({{(Subscript[k, 2,1])}, {(Subscript[k, 2,2])}, {(Subscript[k, 2,3])}, {(Subscript[k, 2,4])}})) = (({{(Subscript[f, 2])((Subscript[x, i]),(Subscript[y, 1,i]), (Subscript[y, 2,i]))}, {(Subscript[f, 2])((Subscript[x, i]) + (1/2)h,(Subscript[y, 1,i]) +  (1/2)(Subscript[k, 1,1]) h , (Subscript[y, 2,i]) +  (1/2)(Subscript[k, 2,1]) h)}, {(Subscript[f, 2])((Subscript[x, i]) +(1/2) h,(Subscript[y, 1,i]) + (1/2)(Subscript[k, 1,2]) h  , (Subscript[y, 2,i]) +  (1/2)(Subscript[k, 2,2]) h)}, {(Subscript[f, 2])((Subscript[x, i]) + h,(Subscript[y, 1,i]) + (Subscript[k, 1,3]) h  , (Subscript[y, 2,i]) + (Subscript[k, 2,3]) h)}}))
Now, we see how convenient it is to deal everything as a vector:
(({{(Subscript[k, 1,1]), (Subscript[k, 2,1])}, {(Subscript[k, 1,2]), (Subscript[k, 2,2])}, {(Subscript[k, 1,3]), (Subscript[k, 2,3])}, {(Subscript[k, 1,4]), (Subscript[k, 2,4])}})) = (({{(Subscript[f, 1])(x; (Subscript[Y, i])), (Subscript[f, 2])(x; (Subscript[Y, i]))}, {(Subscript[f, 1])(x+(1/2) h; (Subscript[Y, i])+(1/2) (Subscript[K, 1]) h), (Subscript[f, 2])(x+(1/2) h; (Subscript[Y, i])+(1/2) (Subscript[K, 1]) h)}, {(Subscript[f, 1])(x+(1/2) h; (Subscript[Y, i])+(1/2) (Subscript[K, 2]) h), (Subscript[f, 2])(x+(1/2) h; (Subscript[Y, i])+(1/2) (Subscript[K, 1]) h)}, {(Subscript[f, 1])(x+ h; (Subscript[Y, i])+ (Subscript[K, 3]) h), (Subscript[f, 2])(x+(1/2) h; (Subscript[Y, i])+(1/2) (Subscript[K, 1]) h)}}))

3.2.1 Simple code

Manipulate[
If[!TrueQ[{lynxold,hareold,hsold}== {lynx,hare,hs}],
(({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}})) =(({{lynx}, {hare}}));
 x = 0. ; np=Ceiling[25./hs];
np=Min[1000,np];h=25/np;
lynxinit={x, (Subscript[y, 1 old])};hareinit={x, (Subscript[y, 2 old])}; 
(Subscript[Y, old])=(({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}})) ;
{lynxPop,harePop}=Reap[For[i=1, i≤ np, i++,{
		(({{(Subscript[k, 11])}, {(Subscript[k, 21])}})) = Flh[x,(Subscript[Y, old])] ;
		(({{(Subscript[k, 12])}, {(Subscript[k, 22])}})) = Flh[x+(h/2),(Subscript[Y, old])+(({{(Subscript[k, 11])}, {(Subscript[k, 21])}})) *(h/2)] ;
         (({{(Subscript[k, 13])}, {(Subscript[k, 23])}})) = Flh[x+(h/2),(Subscript[Y, old])+(({{(Subscript[k, 12])}, {(Subscript[k, 22])}})) *(h/2)] ;
		 (({{(Subscript[k, 14])}, {(Subscript[k, 24])}})) = Flh[x+h,(Subscript[Y, old])+(({{(Subscript[k, 13])}, {(Subscript[k, 23])}})) *h] ;
	(({{(Subscript[y, 1 new])}, {(Subscript[y, 2 new])}})) =  (({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}})) + (h/6) *  (({{(Subscript[k, 11])+2(Subscript[k, 12])+2 (Subscript[k, 13])  + (Subscript[k, 14   ])}, {(Subscript[k, 21])+2(Subscript[k, 22])+2 (Subscript[k, 23])  + (Subscript[k, 24])}})) ;
	x  = x +h ;
		Sow[{ x ,(Subscript[y, 1 new]) },A] ;
		Sow[{ x ,(Subscript[y, 2 new]) },B] ;
		    (({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}}))  =(({{(Subscript[y, 1 new])}, {(Subscript[y, 2 new])}})) ;(Subscript[Y, old]) =  (({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}}))  
			 }]][[2]];
PrependTo[lynxPop,lynxinit];PrependTo[harePop,hareinit];
maxpop=Max[{Transpose[lynxPop][[2]], Transpose[harePop][[2]]}];
LHplot=ListPlot[{lynxPop,harePop}, Joined->True,
PlotStyle->{{Hue[0.,1.,1,0.75], Thickness[0.008]},
		                     {Hue[0.5,1,0.5,0.8], Thickness[0.008]}}
,PlotRange-> {0,maxpop}
,ImageSize->600, AspectRatio-> 0.3] ; 
(* {t,xarr}=Transpose[theta]; *)
{lynxold,hareold,hsold}={lynx,hare,hs}
,];

LHplot

,{{lynx,1},0.5,2},{{hare,5},2,10},{{hs,0.1},0.01,2}
,Initialization:>{
Flh[x_,Y_]:=(({{(-0.9*Y[[1]]+0.5*Y[[1]]*Y[[2]])[[1]]}, {(1.3*Y[[2]]-0.7*Y[[1]]*Y[[2]])[[1]]}})) ;
;}]

3.2.2 Output demo

This is the famous lynx/hare population behavior: we start with a high population of hares (teal) which is an abundant source of food for the lynx (red) with a rapid population growth. But as the lynx population grows large, they deplete the hare population, which declines precipitously, leading to lynx starvation. The lynx population declines, giving the hares a chance to rebuild their population and the cycle repeats itself.

3.3 Using NDSolve

In general, we don’t have to write own codes of common algorithm such as R-K because numerous software packgages have available to use. Below is an example with NDSolve:

Manipulate[
popParm={{dlyn,r},{s,dhar}};
popinit={lynx,hare};
If[!TrueQ[{popParmold,popold}=={popParm,popinit}]
,
solLH=NDSolve[{{yL'[t],yH'[t]}==fLH[t,{yL[t],yH[t]},popParm],{yL[0.],yH[0.]}== popinit},{yL,yH},{t,0,30}];
lhPlot=Plot[Release[{yL[t],yH[t]}/. solLH],{t,0.,30}
,PlotRange-> All
,PlotStyle-> {{Hue[0.,1.,1,0.75], Thickness[0.008]},{Hue[0.5,1,0.5,0.8], Thickness[0.008]}}
,ImageSize->600, AspectRatio-> 0.3] ; 
{popParmold,popold}={popParm,popinit}
,];
lhPlot
,{{lynx,1},0.5,2},{{hare,5},2,10}
,{{dlyn,-0.9},-2,-0.1},{{dhar,2.},0.1,5}
,{{r,0.3},0.001,0.8},{{s,-0.8},-1,0.001}
,Initialization:>{
fLH[x_,Y_,parm_]:={parm[[1,1]]*Y[[1]]+parm[[1,2]]*Y[[1]]*Y[[2]]
,parm[[2,2]]*Y[[2]]+parm[[2,1]]*Y[[1]]*Y[[2]]} ;
}]

4. Higher order differential equation

4.1 Basic concept

We see that solving first-order DE (involving only the first-order derivative) by numerical method is straight forward. So if we encounter a higher order equation, like this:
f(x, y,...) (((d^n) y)/(d (x^n))) + g(x, y,..)(((d^(n-1)) y)/(d (x^(n-1)))) + ...  =0
we should convert if possible to first order.
Example, let's look at the Galileo's pendulum problem:
(((d^2) θ)/(d (t^2)))+(g/l) sin(θ)=0: 
which is second order. So we define an intermediate function:
Ω(t) = ((d θ)/(d t)); and now we have a systems of two equations:
((d Ω)/(d t))= - (g/l) sin(θ) ;
 ((d θ)/(d t)) =Ω(t)

4.2 Example: R-K method to solve the pendulum problem

We can use the same algorithm developed for the lynx-hare problem above for this case. We just need to make a different function

4.2.1 Pendulum simple code

Manipulate[
If[!TrueQ[{θold,hsold}=={θinit,hs}],
 (({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}}))  = (({{0.}, {θinit}})); x = 0. ; np=Ceiling[25./hs];
np=Min[1000,np];h=25/np;
omegainit={x, (Subscript[y, 1 old])};thetainit={x, (Subscript[y, 2 old])}; 
(Subscript[Y, old])=(({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}})) ;
{omega,theta}=Reap[For[i=1, i≤ np, i++,{
		(({{(Subscript[k, 11])}, {(Subscript[k, 21])}})) = F[x,(Subscript[Y, old])] ;
		(({{(Subscript[k, 12])}, {(Subscript[k, 22])}})) = F[x+(h/2),(Subscript[Y, old])+(({{(Subscript[k, 11])}, {(Subscript[k, 21])}})) *(h/2)] ;
         (({{(Subscript[k, 13])}, {(Subscript[k, 23])}})) = F[x+(h/2),(Subscript[Y, old])+(({{(Subscript[k, 12])}, {(Subscript[k, 22])}})) *(h/2)] ;
		 (({{(Subscript[k, 14])}, {(Subscript[k, 24])}})) = F[x+h,(Subscript[Y, old])+(({{(Subscript[k, 13])}, {(Subscript[k, 23])}})) *h] ;
	(({{(Subscript[y, 1 new])}, {(Subscript[y, 2 new])}})) =  (({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}})) + (h/6) *  (({{(Subscript[k, 11])+2(Subscript[k, 12])+2 (Subscript[k, 13])  + (Subscript[k, 14   ])}, {(Subscript[k, 21])+2(Subscript[k, 22])+2 (Subscript[k, 23])  + (Subscript[k, 24])}})) ;
	x  = x +h ;
		Sow[{ x ,(Subscript[y, 1 new]) },A] ;
		Sow[{ x ,(Subscript[y, 2 new]) },B] ;
		    (({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}}))  =	(({{(Subscript[y, 1 new])}, {(Subscript[y, 2 new])}})) ;(Subscript[Y, old]) =  (({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}}))  
			 }]][[2]];
PrependTo[omega,omegainit];PrependTo[theta,thetainit];
pendgr=ListPlot[{omega,theta}, Joined->True,
PlotStyle->{{Hue[0.,1,1], Thickness[0.007]},
		                     {Hue[0.58,1,1], Thickness[0.007]}}
,ImageSize->400] ;
{t,xarr}=Transpose[theta];
{θold,hsold}={θinit,hs}
,];
{pendgr,
Manipulate[
Graphics[{Hue[0.5,1.,1.],Thickness[0.01],Line[{{0.,0.},{Sin[xarr[[i]] ],-Cos[xarr[[i]]  ]}}]
,White,PointSize[0.02],Point[{0.,0.}],
Red,Disk[{Sin[xarr[[i]] ],-Cos[xarr[[i]]  ]},0.1]
,Inset[NumberForm[(i-1.)*h,2],ImageScaled[{0.9,0.9}]
,BaseStyle-> {24, Yellow}]}
,PlotRange-> {{-1.25,1.25},{-1.25,1.25}}
,Background-> Black]
,{i,1,Length[xarr],1,ControlType-> Animator, DefaultDuration->8
,AnimationRunning-> False}]
}

,{{θinit,Pi/2},0.01,0.9 Pi},{{hs,0.01},0.005,0.5}
,Initialization:>{
F[t_,Y_List]:= (({{-Sin[Y[[2]] ][[1]]}, {Y[[1]][[1]]}}));}]

4.2.2 output

4.3 Pendulum with NDSolve

4.4 Example: R-K method for simple RLC circuit (damped oscillation)

Below is a circuit

With Laplace transform, we can determine the behavior based on the poles of the transfer function.

Underdamped, complex roots near Im axis, highly oscillatory response

Overdamped, real negative roots,  rapid exponential decay response

Of course we know the solution analytically, but we use this example to show how to apply Runge-Kutta and see how it works.

Let the general equation be:
 (((d^2) y)/(d (t^2)))+2 ζ (Subscript[ω, n])((dy)/(dt))+((Subscript[ω, n])^2)y=0: 
which is second order. So we define an intermediate function:
 (Subscript[y, 1])(t)=((dy)/(dt)); and now we have a systems of two equations:
((d(Subscript[y, 1]))/(dt))=-2 ζ (Subscript[ω, n])(Subscript[y, 1])-((Subscript[ω, n])^2)(Subscript[y , 0])
 ((d(Subscript[y, 0]))/(dt)) =(Subscript[y, 1])(t) 
We can use the same algorithm developed for the lynx-hare problem above for this case. We just need to make a different vector function

Manipulate[
If[!TrueQ[{ωnold,ζold,xendold,hsold}=={ωn,ζ,xend,hs}],
h=Max[hs,0.001 xend] ;
np=Ceiling[xend/h];
If[np≤ 4, h=xend/4.;np=Ceiling[xend/h];,];
(({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}}))  = (({{0.}, {1}})); x = 0. ; 
volt0init={x, (Subscript[y, 1 old])};
volt1init={x, (Subscript[y, 2 old])};
(Subscript[Y, old])=(({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}})) ;
{volt0,volt1}=Reap[For[i=1, i≤ np,i++ ,{
(({{(Subscript[k, 11])}, {(Subscript[k, 21])}})) = Fosc[(Subscript[Y, old]),{ωn,ζ}] ;
(({{(Subscript[k, 12])}, {(Subscript[k, 22])}})) = Fosc[(Subscript[Y, old])+(({{(Subscript[k, 11])}, {(Subscript[k, 21])}})) *(h/2),{ωn,ζ}] ;(({{(Subscript[k, 13])}, {(Subscript[k, 23])}})) =Fosc[(Subscript[Y, old])+(({{(Subscript[k, 12])}, {(Subscript[k, 22])}})) *(h/2),{ωn,ζ}] ;(({{(Subscript[k, 14])}, {(Subscript[k, 24])}})) =Fosc[(Subscript[Y, old])+(({{(Subscript[k, 13])}, {(Subscript[k, 23])}})) *h,{ωn,ζ}] ;(({{(Subscript[y, 1 new])}, {(Subscript[y, 2 new])}})) =  (({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}})) + (h/6) *  (({{(Subscript[k, 11])+2(Subscript[k, 12])+2 (Subscript[k, 13])  + (Subscript[k, 14   ])}, {(Subscript[k, 21])+2(Subscript[k, 22])+2 (Subscript[k, 23])  + (Subscript[k, 24])}})) ;
x  = x +h ;
Sow[{x ,(Subscript[y, 1 new])},V1];Sow[{x ,(Subscript[y, 2 new])},V2] ; 
(({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}}))  =(({{(Subscript[y, 1 new])}, {(Subscript[y, 2 new])}})) ;(Subscript[Y, old]) =(({{(Subscript[y, 1 old])}, {(Subscript[y, 2 old])}}))  
			 }]][[2]];
PrependTo[volt0,volt0init];PrependTo[volt1,volt1init];

{ωnold,ζold,xendold,hsold}={ωn,ζ,xend,hs}
,];
ListPlot[{volt0,volt1}, Joined-> True
,PlotStyle->{{RGBColor[1.,0.,0], Thickness[0.007]},
		                     {RGBColor[0.,0.,1], Thickness[0.007]}}
,Frame-> True, GridLines-> Automatic
,FrameLabel-> {"time","voltage"}
,LabelStyle-> {16,Bold,Black,FontFamily-> "Arial"}
,ImageSize-> {400.,300}
,PlotRange->{-1.1,1.1}]
,Grid[{{Control[{{ωn,1.,Button["nat. frequency",BaseStyle->{14,Bold, FontFamily-> "Calibri"},Background->Hue[0.8,0.4,1]]},0.5,5,0.1}]
,Control[{{xend,20,Button["end time",BaseStyle->{14,Bold, FontFamily-> "Calibri"},Background->Hue[0.25,0.8,1]]},1,40,0.025}]}
,{Control[{{ζ,0.2,Button["damping coef.",BaseStyle->{14,Bold, FontFamily-> "Calibri"},Background->Hue[0.07,0.8,1]]},0.001,2.,0.025}]
,Control[{{hs,0.1,Button["step size",BaseStyle->{14,Bold, FontFamily-> "Calibri"},Background->Hue[0.5,0.4,1]]},0.01,0.5}]
}}]
,Initialization:>{
Fosc[Y_,{ωn_,ζ_}]:=(({{Y[[2,1]]}, {-ωn (ωn Y[[1,1]]+2ζ Y[[2,1]])}})) ;

}]