x = seq(-2*pi, 2*pi, length.out = 250)
y = sin(x)

plot(x,y,
     type = "l",
     col = "grey",
     lwd = 2,
     bty = "n")
x0 = pi/2

pol_sin = function(x, x0, n){
  y = sin(x0);
  for(i in 1:n){
    if(i%%4 == 0){
      y = y + sin(x0) * (x-x0)^i/factorial(i)
    } else if(i%%4 == 1) {
      y = y + cos(x0) * (x-x0)^i/factorial(i)
    } else if(i%%4 == 2) {
      y = y - sin(x0) * (x-x0)^i/factorial(i)

    } else if(i%%4 == 3) {
      y = y - cos(x0) * (x-x0)^i/factorial(i)
    }
  }
  return(y)
}
pol_sin = Vectorize(pol_sin, vectorize.args = "x")
lines(x,pol_sin(x,x0,2),
      col = "red",
      lwd = 2, lty = 2)
lines(x,pol_sin(x,x0,3),
      col = "blue",
      lwd = 2, lty = 2)
lines(x,pol_sin(x,x0,6),
      col = "pink",
      lwd = 2, lty = 2)
lines(x,pol_sin(x,x0,12),
      col = "magenta",
      lwd = 2, lty = 2)
lines(x,pol_sin(x,x0,20),
      col = "forestgreen",
      lwd = 2, lty = 2)