dat<-data.frame(x=c(2,3,4,5,6),y=c(4,4,6,6,10))

N = length(dat$x)
x = dat$x-mean(dat$x)
y = dat$y-mean(dat$y)

library(bbmle)
fn<-function(b, sigma) {
  (N/2)*log(sigma^2)+1/(2*sigma^2)*(sum((y-b*x)^2))
}
mle<-mle2(fn,start=list(b=1.3, sigma=1))
summary(mle)

Lik <- function(param1, param2, x, y, N){
  beta1 = param1  
  sigma = param2  
  Likelihood = vector("numeric",length(beta1))
  for (i in 1:length(beta1)){
    Likelihood[i] = (2*pi*sigma[i]^2)^(-N/2)*exp((-1/(2*sigma[i]^2))*sum((y-beta1[i]*x)^2))
  }
  return(Likelihood)
}

param1 <- seq(1.1, 2.1, by=0.05)
names(param1) <- param1
param2 <- seq(0.5, 1.0, by=0.02)
names(param2) <- param2

z <- outer(X = param1, Y = param2, FUN = Lik, N = N, x = x, y = y)
max(z)
zd<-data.frame(z)

inds = which(zd == max(zd), arr.ind=TRUE)
inds
(rnames = rownames(zd)[inds[,1]])
(cnames = colnames(zd)[inds[,2]])

require(lattice)
#> Loading required package: lattice
wireframe(z, drape=T, xlab="beta", ylab="sigma", zlab="Lik",col.regions=rainbow(100))
wireframe(z, xlab="beta1", ylab="sigma", zlab="Lik")

# test it out
Lik(param1[3], param2[3],  x = x, y = y, N = N)
z[3,3]