if( !require('phaseR') ){ install.packages('phaseR') } library(phaseR) my_yprime = function(t, y, parameters){ a = parameters[1] b = parameters[2] y_prime = rep( NA, length(y) ) lt_zero_mask = y <= 0 y_prime[lt_zero_mask] = 0 gt_zero_mask = !lt_zero_mask y_prime[gt_zero_mask] = a*y[gt_zero_mask] - b*sqrt(y[gt_zero_mask]) list( y_prime ) } diff_eq_params = c( 2, 4 ) #postscript("../../WriteUp/Graphics/Chapter2/chap_2_sect_5_prob_11_plot.eps", onefile=FALSE, horizontal=FALSE) par(mfrow=c(1, 2)) ys = seq( 0.1, 5, length.out=100 ) fys = my_yprime( 0.0, ys, diff_eq_params ) fys = fys[[1]] plot( ys, fys, type='l', xlab='y', ylab='f(y)' ) abline(h=0, col='black') grid() # Plot the flow field, nullclines and several trajectories: # flowField(my_yprime, x.lim = c(0, 4), y.lim = c(0, 6), parameters = diff_eq_params, points = 21, system = "one.dim", add = FALSE, xlab='t') abline(h=c(0, 4), col='red') trajectory(my_yprime, y0 = c(3.5, 4.5), t.end = 4, parameters = diff_eq_params, system = "one.dim") par(mfrow=c(1, 1)) #dev.off()