Can't reproduce confint profile likelihood results

Question

I'm trying to follow this answer here and the linked paper. But I can't seem to get the same results as the confint gives. Sometimes I get something very simmilar, but sometimes my CI is more conservative or less conservative than confint.

I want to know also if the general idea of the implementation is correct - say your model is $a+ b*x$, and you wish to have a CI for $b$ - are you supposed to find, for every value of $b$ the maximum log-likelihood of $a$, and compute these log-likelihoods, and find the $b$'s whose log-likelihood is bigger than $ll(b^{MLE}) - \chi^2_{1,0.975}/2$ ?

Here's the code:

# Generate the data
rpt = 4
x = rep(seq(1, 10, 0.5), times=rpt)
mu = 2.3*x # + rnorm(length(x), sd=0.5)
y = rpois(length(mu), mu)
plot(x, y)
fit = glm(y ~ x, family=poisson(link="identity"))
fit$coefficients
profile log likelihood function
loglik = function(b) {
  func = function(a) {
    mu.hat = a + b*x
    ll = 0
    for (i in 1:length(x)) {
      ll = ll + dpois(y[i], mu.hat[i], log=T)
    }
    return(ll)
  }
  fit1 = optimize(func, c(0, 5), maximum=T)
  return(fit1$objective)
}
logLik(fit) # built-in function that computes the log-likelihood
loglik(fit$coefficients[2]) # profile ll - should be equal for MLE
beta = seq(1.9, 2.7, 0.0001)
ll = rep(0, length(beta))
for (i in 1:length(beta)){
  ll[i] = loglik(beta[i])
}
plot(beta, ll)
wch = which(ll > logLik(fit)-qchisq(0.975, 1)/2)
(ind.lo = wch[1])
(beta.lo = beta[ind.lo])
(ind.hi = wch[length(wch)])
(beta.hi = beta[ind.hi])
confint(fit) # compare to confint
(beta.MLE = fit$coefficients[2])
pchisq(2(loglik(beta.MLE)-loglik(beta.lo)), 1)
pchisq(2(loglik(beta.MLE)-loglik(beta.hi)), 1)

UPDATE:

Using this glm seems to give exact results:

loglik = function(b) {
  fit1 = glm(y ~ 1, offset=b*x, family=poisson(link="identity"))
  return(logLik(fit1))
}