#postscript("../../WriteUp/Graphics/Chapter2/chap_2_sect_3_prob_15_plot.eps", onefile=FALSE, horizontal=FALSE) k = 1/5 ts = seq( 0, 20, length.out=100 ) # Compute the numerical integral part: # library(stats) integrand_f = function(t){ exp( -t/5 + cos(t)/5 ) } int_part = rep( NA, length(ts) ) for( ti in 1:length(ts) ){ res = integrate( integrand_f, 0, ts[ti] ) int_part[ti] = res\$value } y0 = 1/2 y_1 = y0 * exp( ts/5 - cos(ts)/5 + 1/5 ) - k * exp( ts/5 - cos(ts)/5 ) * int_part y0 = 3/4 y_2 = y0 * exp( ts/5 - cos(ts)/5 + 1/5 ) - k * exp( ts/5 - cos(ts)/5 ) * int_part y0 = 1 y_3 = y0 * exp( ts/5 - cos(ts)/5 + 1/5 ) - k * exp( ts/5 - cos(ts)/5 ) * int_part y_min = min( c( y_1, y_2, y_3 ) ) y_max = max( c( y_1, y_2, y_3 ) ) plot( ts, y_1, ylim=c( y_min, y_max ), col='blue', type='l', xlab='t', ylab='y(t)' ) lines( ts, y_2, col='green', type='l' ) lines( ts, y_3, col='red', type='l' ) abline(h=0, col='black') legend( 'bottomleft', c('y0=1/2', 'y0=3/4', 'y0=1'), lty=1, col=c('blue', 'green', 'red') ) grid() #dev.off() # Part (b): # rhs = k*exp(-1/5)*int_part plot( ts, rhs, type='l', xlab='t', ylab='rhs', main='integral part' ) print( rhs[length(rhs)] )