fx <- expression( sin(x) ); f <- function(x) eval(fx)
g1x <-D(fx,"x"); g1 <- function(x) eval(g1x)
g2x <-D(g1x,"x"); g2 <- function(x) eval(g2x)
g3x <-D(g2x,"x"); g3 <- function(x) eval(g3x)
g4x <-D(g3x,"x"); g4 <- function(x) eval(g4x)
g5x <-D(g4x,"x"); g5 <- function(x) eval(g5x)
g6x <-D(g5x,"x"); g6 <- function(x) eval(g6x)
g7x <-D(g6x,"x"); g7 <- function(x) eval(g7x)
g8x <-D(g7x,"x"); g8 <- function(x) eval(g8x)
g9x <-D(g8x,"x"); g9 <- function(x) eval(g9x)
g10x <-D(g9x,"x"); g10 <- function(x) eval(g10x)
g11x <-D(g10x,"x"); g11 <- function(x) eval(g11x)
g12x <-D(g11x,"x"); g12 <- function(x) eval(g12x)
g13x <-D(g12x,"x"); g13 <- function(x) eval(g13x)
p1 <- function(x) f(0)+g1(0)*x
p2 <- function(x) p1(x)+g2(0)*x^2/factorial(2)
p3 <- function(x) p2(x)+g3(0)*x^3/factorial(3)
p4 <- function(x) p3(x)+g4(0)*x^4/factorial(4)
p5 <- function(x) p4(x)+g5(0)*x^5/factorial(5)
p6 <- function(x) p5(x)+g6(0)*x^6/factorial(6)
p7 <- function(x) p6(x)+g7(0)*x^7/factorial(7)
p8 <- function(x) p7(x)+g8(0)*x^8/factorial(8)
p9 <- function(x) p8(x)+g9(0)*x^9/factorial(9)
p10 <- function(x) p9(x)+g10(0)*x^10/factorial(10)
p11 <- function(x) p10(x)+g11(0)*x^11/factorial(11)
p12 <- function(x) p11(x)+g12(0)*x^12/factorial(12)
p13 <- function(x) p12(x)+g13(0)*x^13/factorial(13)
c <- c(1,4,3,2,6,4,3,2,6,4,3,2,6); j <- 1
a <- -3; b <- 7
plot(f,a,b,ylim=c(-1.7,1.7), col=c[j])
abline(h=seq(-1.5,1.5,0.5),v=seq(a,b,1),lty=3); abline(v=0,h=0,lty=2)
points(0,0)
j <- j+1; plot(p1,a,b,add=TRUE,col=c[j])
j <- j+1; plot(p3,a,b,add=TRUE,col=c[j])
j <- j+1; plot(p5,a,b,add=TRUE,col=c[j])
j <- j+1; plot(p7,a,b,add=TRUE,col=c[j])
j <- j+1; plot(p9,a,b,add=TRUE,col=c[j])
j <- j+1; plot(p11,a,b,add=TRUE,col=c[j])
j <- j+1; plot(p13,a,b,add=TRUE,col=c[j])