a1 <- c(1, 1, 1)
a2 <- c(2, 2, 2)
b <- c(-1, 0, 1)
theta <- seq(-3, 3, 0.1)

# calculate P for each item 
p1 <- matrix(NA, length(theta), length(b)) # create a matrix of P
for (i in 1:length(b)){
  p1[,i] <- 1 / (1 + exp(- a1[i] * (theta - b[i])))    
}

p2 <- matrix(NA, length(theta), length(b)) # create a matrix of P
for (i in 1:length(b)){
  p2[,i] <- 1 / (1 + exp(- a2[i] * (theta - b[i])))    
}

lik1 <- p1[,1]*p1[,2]*(1-p1[,3]) # 110
lik2 <- p2[,1]*p2[,2]*(1-p2[,3])

plot(theta, lik1, type="l", ylim=range(lik1, lik2),
     xlab=expression(paste("Latent trait (",theta, ")")), 
     ylab="Likelihood")
lines(theta, lik2, lty=2)
legend("topleft", lty=1:2, legend = c("a=1", "a=2"))