# X0     : Punto inicial 0
# X1     : Punto inicial 1
# f(x)   : La función evaluada en el punto x

# n      : La iteración actual
# It     : Máximo numero de Iteraciones

# Fórmula Utilizada
#  m = (f(X[n+1]) - f(X[n]))/(X[n+1] - X[n])
#  X[n+2] = X[n+1] - f(x[n+1])/m
#La ventaja sobre secante es que el nuevo punto cae entre los 2 anteriores

#La tabla tiene las siguientes columnas
#iteración - x - f(x) - (x1-x0) - [f(x1)-f(x0)] - 1/m - f(x)/m

# esta es la función a la que deseamos calcularle la raíz
f = function(x) {
  3*x^3-4*x^2+5*x-6
}

#secante utiliza 2 puntos iniciales.
#Realizaremos el método usando la tabla.
posifalsa = function(x0=3,x1=2,It=100) {

  tabla = matrix(c(1,x0,f(x0),0,0,0,0),nrow=1)    #La primera línea de la tabla con el punto x0
  tabla = rbind(tabla,c(2,x1,f(x1),0,0,0,0))      #Agregamos la segunda línea con X1
  tabla[2,4] = x1 - x0                            #La diferencia entre los puntos
  tabla[2,5] = f(x1) - f(x0)                      #La diferencia entre las funciones
  tabla[2,6] = tabla[2,4]/tabla[2,5]              #Calculamos 1/m con la formula punto pendiente
  tabla[2,7] = tabla[2,6]*f(x1)                   #Similar a Newton Hallamos f(x0)/m

  n = 2                                           #Tenemos 2 filas de la tabla hechas hasta el momento
  x = x1                                          #El valor del último x
  p = n - 1                                       
  while(abs(f(x))>10^-4 && n<=It) {               #Condición de parar o seguir
    n = n + 1                                     #pasamos a la línea siguiente
    x = tabla[p,2] - tabla[n - 1,7]               #Calculamos el nuevo x con la formula x = x - f(x)/m
    if (f(x)*tabla[p,3]<0) 
       p = n-1   
    else
       p = n-2
    
    tabla = rbind(tabla,c(n,x,f(x),0,0,0,0))      #agregamos un línea a la tabla con el nuevo x
    tabla[n,4] = tabla[n,2] - tabla[n-1,2]        #Calculamos diferencias de x y el anterior
    tabla[n,5] = tabla[n,3] - tabla[n-1,3]        #Calculamos diferencias de f(x) y f anterior
    tabla[n,6] = tabla[n,4] / tabla[n,5]          #Dividimos para calcular 1/m
    tabla[n,7] = tabla[n,6] * tabla[p,3]          #multiplicamos para hallar f(x)/m
  }#while
  tabla = as.data.frame(tabla)                    #convertimos la matriz en datos.
  names(tabla) = c("iteracion","x1","f(x1)","x1-X0","f(x1)-f(X0)","1/m","f(x1)/m") # Le ponemos nombre a las columnas
  tabla                                           #Mostramos o retornamos la tabla 
}

posifalsa(3,1)                                      #Probamos la función