# Solve a set of second-order differential equations

#dy/dt = F(x,y,t) y=ω
#dx/dt =G(x,y,t)　x=θ

# gravitational acceleration g [m/s^2]
g=9.8

# rod length　L [m]
L=1.000

# mass of the weight m [kg]
m=0.15

# viscous damping coefficient c [Ns]
c=0.000464

# right side of differential equation
F(x,y,t)= -(g/L)*sin(x)-c/(m*L)*y
G(x,y,t)=y

# time width and the step numbers of Runge-Kutta step
dt = 0.001
N = 10000

# default value of x, x0 [rad] 
x0=0.03

# default value of y=dx/dt y0 [rad/s]
y0=0

# initializing the variable of Runge-Kutta method
k1=0.0; k2=0.0; k3=0.0; k4=0.0
m1=0.0; m2=0.0; m3=0.0; m4=0.0

outputfile="differential-eq-output.dat"

set samples N
set xrange [0:dt*(N-1)]

# writing down the data to separate file
set table outputfile
plot xold=x0,yold=y0,\
"+" using ($1+dt):\
(k1=dt*F(xold,yold,$1), \
m1 =dt*G(xold,yold,$1), \
k2 =dt*F(xold+m1*0.5,yold+k1*0.5, $1+dt*0.5), \
m2 =dt*G(xold+m1*0.5,yold+k1*0.5, $1+dt*0.5), \
k3 =dt*F(xold+m2*0.5,yold+k2*0.5, $1+dt*0.5), \
m3 =dt*G(xold+m2*0.5,yold+k2*0.5, $1+dt*0.5), \
k4 =dt*F(xold+m3,yold+k3, $1+dt), \
m4 =dt*G(xold+m3,yold+k3, $1+dt), \
ynew=yold+(1/6.0)*(k1+2.0*k2+2.0*k3+k4),\
xnew=xold+(1/6.0)*(m1+2.0*m2+2.0*m3+m4),\
xold=xnew, \
yold=ynew, \
xnew)

unset table

# designating the range of x and y axes
set xrange [0:10]
set yrange [-0.05:0.05]

#  setting the label of x and y axes
set xlabel "t [sec]"
set ylabel "θ [rad]"

# graphing
plot outputfile u 1:2 w l linecolor -1 linewidth 2 \
notitle,0 w l linecolor -1 notitle

# designating the output format and option
set terminal jpeg

# designating the output file name
set output 'rk4.jpeg'

# regraphing
replot

# resetting the format and output
set output
set terminal win
[EOF]