# Real dumped pendulum: small oscillation (pend_re_soN).
# as Code 3.1 Ideal pendulum: small oscillations (pend_id_s_o.R),
# but here with the dumping constant gamma=3 (towards a node);
# see also Real dumped pendulum: small oscillation (pend_re_soF) (focus).

parms<- c(1,2,3)
tinit<-0
tfin<-10
step<-0.01
times<-seq(tinit,tfin,by=step)
funct<- function(t,integ,p){ # the system of equations
x<-integ[1]
y<-integ[2]
dx<-  parms[1]*y
dy<- -parms[3]*y-parms[2]*x
list(c(dx,dy))
} # end of funct
require(deSolve)
cinit<- c(2.5,2.5)
xy<-lsoda(cinit,times,funct,parms)
# to plot different trajectories, you have to choose:
# either in phase space or in time domain. With two runs you have both.

# phase space:
plot(xy[,2],xy[,3],type="l",xlab="x(t)",ylab="y(t)",xlim=c(-3,3),ylim=c(-3,3),
cex.lab=1.5,cex.axis=1.2,lwd=3,lty=1)
# time domain:
#plot(xy[,1],xy[,2],type="l",xlab="t",ylab="x(t)",xlim=c(0,10),ylim=c(-3,3),
#cex.lab=1.5,cex.axis=1.2,lwd=3,lty=1)
## second  initial point
cinit<-c(-2.5,-2.5)
xy<- lsoda(cinit,times,funct,parms)
lines(xy[,2],xy[,3],col="red",lwd=3,lty=2)   # phase space
#lines(xy[,1],xy[,2],col="red",lwd=3,lty=2)  # time domain

# further initial points
cinit<-c(1.5,1.5)
xy<- lsoda(cinit,times,funct,parms)
lines(xy[,2],xy[,3],col="magenta",lwd=4,lty=3)   # phase space
#lines(xy[,1],xy[,2],col="magenta",lwd=8,lty=3)  # time domain

cinit<-c(-1.5,-1.5)
xy<- lsoda(cinit,times,funct,parms)
lines(xy[,2],xy[,3],col="blue",lwd=3,lty=4)   # phase space
#lines(xy[,1],xy[,2],col="blue",lwd=3,lty=4)  # time domain