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

parms<- c(1,2,0.25)
tinit<-0
tfin<-50
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,3.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(-4,4),ylim=c(-4,4),
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,50),ylim=c(-3,3),
#cex.lab=1.5,cex.axis=1.2,lwd=3,lty=1)
## second  initial point
cinit<-c(2,0.2)
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