Wagner, veja a adaptação de seu código comentada abaixo.
# função
Lmax_fun <- function(par){
L <- par
Ho <- 12
Hvar <- 0.09061577
So <- -0.03
k <- 4.182488e-06
((Ho*Hvar/((-So-(k*L^1.75))*(1-Hvar)))-L)^2
}
L <- 1:300
Ho <- 12; Hvar <- 0.09061577; So <- -0.03; k <- 4.182488e-06
DEN_L <- ((-So-(k*L^1.75))*(1-Hvar))*L
plot(L,DEN_L,type = "l", ylim = c(0,max(DEN_L)*1.10))
abline(h=Ho*Hvar)
# otimização unidimensional - optimize()
# resultado depende do intervalo escolhido
# pelo gráfico da função podemos ver os intervalos das raízes
optimize(f = Lmax_fun, interval = c(min(L), max(L)))
optimize(f = Lmax_fun, interval = c(1, 200))
optimize(f = Lmax_fun, interval = c(100, 300))
library(pracma)
# findmins() divide o intervalo de busca n vezes (n = 100 por default)
mins <- findmins(f = Lmax_fun, min(L), max(L))
points(mins, rep(Ho*Hvar, 3), pch = 20)
#Condição para valer a escolha da função
mins[abs(So)/(k*mins^1.75) >= 1.75+1]