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`

<!-- language: r -->

    summary.lm
    #R ...
    #R    rdf <- z$df.residual 
    #R    ... 
    #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 <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
    #R    se <- sqrt(diag(R) * resvar)
    #R    ...

This is the [observed Fisher information](https://en.wikipedia.org/wiki/Observed_information) for $(\beta_0, \beta_1)$ conditional on $\hat\sigma^2$. Now the 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

<!-- language: r -->

    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