one.r

## a second order differential equation
## D2 y + 2 r w D y + w^2 y = 0
# D is the differential operator
# r is the damping ratio (r = c/(2mw))
# w is the undamped frequency
# c is the viscous damping coefficent
# m is the mass
#////////////////////////////////////////////////////////////////////////////////
#//                                                                            //
#//  Description:                                                              //
#//     A differential equation of order greater than one can be transformed   //
#//     to a multidimensional first order differential equation.  In           //
#//     particular, given a second order differential equation y'' = f(x,y,y') // 
#//     define z = y', then the pair (y,z)' = (z, f(x,y,z) ) is a first order  //
#//     differential equation.  The initial conditions y(x0) = y[0],           //
#//     y'(x0) = c becomes the initial condition (y,z)(x0) = (y[0],c).         //
#//     when x = x0 numerically evaluates f(x,y) four times per step. Then for //
#//     step i+1,                                                              //
#//           y[i+1] = y[i] + h * ( yp[i] + 1/6 (k1 + k2 + k3 ) and            //
#//           yp[i+1] = yp[i] + 1/6 * ( k1 + 2 k2 + 2 k3 + k4 ) where          //
#//         k1 = h * f(x[i], y[i], yp[i]),                                     //
#//         k2 = h * f(x[i]+h/2, y[i]+(h/2)*yp[i], yp[i]+k1/2),                //
#//         k3 = h * f(x[i]+h/2, y[i]+(h/2)*yp[i]+(h/4)*k1, yp[i]+k2/2),       //
#//     x[i] = x0 + i * h.                                                     //
#//                                                                            //
#////////////////////////////////////////////////////////////////////////////////

f <- function(t_n,y_n,yp_n, r, w) {
# ## D2 y + 2 r w D y + w^2 y = 0 or # f = f(t, y, yp )
	ans = -2.0*r*w*yp_n - w*w*y_n
	ans
}

### Runge Kutta 2order, call every step
## incomplete!
## incomplete!
## incomplete!

integrate_f <- function(t_n,y_n,yp_n,h=0.0009765, r=1.0, w = 1.0) {
	one_sixth = 1.0/6.0
	half_h = 0.5 *h
	
	k1 = h*f(t_n, y_n, yp_n, r, w)
	k2 = h*f(t_n + half_h, y_n + half_h*y_p, yp_n + 0.5*k1, r, w)
	k3 = h*f(t_n + half_h, y_n + half_h*(y_p + 0.5*half_h*k1), yp_n + 0.5*k2, r, w)
	k4 = h* f(t_n + h, y_n + yp_n*h + half_h*k2, yp_n + k3)
	
	y_n_new = y_p + h*(yp_n + one_sixth*(k1+k2+k3))
	yp_n_new = yp_n + one_sixth*(k1+2*k2+2*k3+k4)
}

two.r

##declare equation F(t,x,v)= -(k/m)v-(b/m)x
input k,b,m 
x(0) <- 0 ##declare initial 
v(0) <- 0 ##declare initial velocity
h <- 0.2 ##declare time step 
ti <- 0 ## declare tinitial
tf <- 5 ## declare tfinal
x <- x(0)
v <- v(0)
loop for t from ti to tf with increment of h

x <-array(, dim = c(25,2))
		dx1 = h*v
		dv1 = h*F(t,x,v)
		
		dx2 = h*(v+dv1/2)
		dv2 = h*F(t+h/2, x+dx1/2, v+dv1/2)
		
		dx3 = h*(v+dv2/2)
		dv3 = h*F(t+h/2, x+dx2/2, v+dv2/2)
		
		dx4 = h*(v+dv3)
		dv4 = h*F(t+h x+dx3, v+dv3)
		
		dx = (dx1 + 2*dx2 + 2*dx3 + dx4)/6
		dv = (dv1 + 2*dv2 + 2*dv3 + dv4)/6
		
		x = x+dx
		v = v+dv
end loop