if( !require('phaseR') ){ install.packages('phaseR') } library(phaseR) A = matrix(c(-1, 0, 0, -1), nrow=2, ncol=2, byrow=T) my_yprime = function(t, y, parameters) { dy = parameters$A %*% y list(dy) } diff_eq_params = list(A=A) #postscript('../../WriteUp/Graphics/Chapter9/chap_9_sect_1_prob_11_plot.eps', onefile=FALSE, horizontal=FALSE) par(mfrow=c(1, 2)) t.end = 3.2 flowField(my_yprime, x.lim = c(-1, +1), y.lim = c(-1, +1), parameters = diff_eq_params, points = 21, add = FALSE, xlab='x_1', ylab='x_2') trajectory(my_yprime, y0 = c(-0.2, 0.5), t.end = t.end, parameters = diff_eq_params, col='black', pch=19) trajectory(my_yprime, y0 = c(0.2, -0.75), t.end = t.end, parameters = diff_eq_params, col='red', pch=19) ts = seq(0, t.end, length.out=100) res_1 = ode(y=c(-0.2, 0.5), times=ts, fun=my_yprime, parms=diff_eq_params) x_1_1 = res_1[,'1'] res_3 = ode(y=c(0.2, -0.75), times=ts, fun=my_yprime, parms=diff_eq_params) x_1_3 = res_3[,'1'] x_min = min( c( x_1_1, x_1_3) ) x_max = max( c( x_1_1, x_1_3) ) plot(ts, x_1_1, ylim=c(x_min, x_max), type='l', xlab='time', ylab='x_1(t)', col='black') points(ts, x_1_3, type='l', col='red') grid() par(mfrow=c(1, 1)) #dev.off()