The issues is that the standard errors comes from
$$\hat\sigma^2 (X^\top X)^{-1}$$
where $\hat\sigma^2$ is the unbaised estimator and not the MLE. See summary.lm
summary.lm
#R function (object, correlation = FALSE, symbolic.cor = FALSE,
#R ...)
#R {
#R z <- object
#R p <- z$rank
#R rdf <- z$df.residual
#R ...
#R r <- z$residuals
#R f <- z$fitted.values
#R w <- z$weights
#R if (is.null(w)) {
#R mss <- if (attr(z$terms, "intercept"))
#R sum((f - mean(f))^2)
#R else sum(f^2)
#R rss <- sum(r^2)
#R }
#R ...
#R resvar <- rss/rdf
#R ...
#R R <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
#R se <- sqrt(diag(R) * resvar)
#R ...
This is the inverse observed Fisher information for $(\beta_0, \beta_1)$ conditional on $\hat\sigma^2$. Now the inverse observed Fisher information you compute is for the triplet $(\beta_0, \beta_1, \sigma)$. Thus, I gather the standard errors should differ by factor $\sqrt{n/(n-3 + 1)}$ or something similar. This seems to be the case
set.seed(1)
n = 4 # very small sample size !
b0 <- 5
b1 <- 2
sigma <- 5
x <- runif(n, 1, 100)
y = b0 + b1*x + rnorm(n, 0, sigma)
negLL <- function(beta, y, x) {
b0 <- beta[1]
b1 <- beta[2]
sigma <- beta[3]
yhat <- b0 + b1*x
return(-sum(dnorm(y, yhat, sigma, log = TRUE)))
}
res <- optim(c(0, 0, 1), negLL, y = y, x = x, hessian=TRUE)
estimates <- res$par # Parameters estimates
(se <- sqrt(diag(solve(res$hessian))))
#R [1] 5.690 0.097 1.653
k <- 3
se * sqrt(n / (n-k+1))
#R [1] 8.047 0.137 2.338