source( 'utils.R')

t = seq( 0, 9, length.out=500 )

k = 1/2
f_onehalf = ( u_sub_c(t, 4-k) - u_sub_c(t, 4+k) ) / (2*k)

k = 1
f_one = ( u_sub_c(t, 4-k) - u_sub_c(t, 4+k) ) / (2*k)

k = 2
f_two = ( u_sub_c(t, 4-k) - u_sub_c(t, 4+k) ) / (2*k)

k = 3
f_three = ( u_sub_c(t, 4-k) - u_sub_c(t, 4+k) ) / (2*k)

#postscript("../../WriteUp/Graphics/Chapter6/sect_4_prob_18_f_k_plots.eps", onefile=FALSE, horizontal=FALSE)
plot( t, f_onehalf, type='l', col='black', xlab='t', ylab='f_k', main='f_k for various k' )
lines( t, f_one, type='l', col='blue' )
lines( t, f_two, type='l', col='red' )
lines( t, f_three, type='l', col='green' )
grid()
#dev.off()


# Plot the solution:
#
t = seq( 0, 40, length.out=1000 )
f = function(t){
    r = sqrt(143)/6
    1/4 - (1/4) * exp( -t/6 ) * cos( r * t ) - (1/(4*sqrt(143))) * exp( -t/6 ) * sin( r * t )
}


k = 1/2
y_onehalf = (1/(2*k)) * ( u_sub_c(t, 4-k ) * f( t-4+k ) - u_sub_c(t, 4 +k ) * f( t-4-k ))
k = 1
y_one = (1/(2*k)) * ( u_sub_c(t, 4-k ) * f( t-4+k ) - u_sub_c(t, 4+k ) * f( t-4-k ) )
k = 2
y_two = (1/(2*k)) * ( u_sub_c(t, 4-k ) * f( t-4+k ) - u_sub_c(t, 4+k ) * f( t-4-k ) )

#postscript("../../WriteUp/Graphics/Chapter6/sect_4_prob_18_yt_plots.eps", onefile=FALSE, horizontal=FALSE)
plot( t, y_onehalf, type='l', col='black', xlab='t', ylab='y(t)' )
lines( t, y_one, type='l', col='blue' )
lines( t, y_two, type='l', col='red' )
legend( 'bottomright', c('k=1/2', 'k=1', 'k=2'), lwd=2, lty=c(1, 1, 1), col=c('black', 'blue', 'red') )
grid()
#dev.off()