I am trying to predict a binary GLM with multiple predictors. I can do it fine with one predictor variable however struggle when I use multiple

Sample data:

structure(list(attempt = structure(c(1L, 2L, 1L, 2L, 1L, 1L, 
1L, 2L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 
2L, 1L, 2L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 
1L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 1L, 2L, 1L, 2L), .Label = c("1", 
"2"), class = "factor"), searchtime = c(137, 90, 164, 32, 39, 
30, 197, 308, 172, 48, 867, 117, 63, 1345, 38, 122, 226, 397, 
0, 106, 259, 220, 170, 102, 46, 327, 8, 10, 23, 108, 315, 318, 
70, 646, 69, 97, 117, 45, 31, 64, 125, 17, 240, 63, 549, 1651, 
233, 406, 334, 168, 127, 47, 881), mean.search.flow = c(15.97766667, 
14.226, 17.15724762, 14.7465, 39.579, 23.355, 110.2926923, 71.95709524, 
72.73666667, 32.37466667, 50.34905172, 27.98471429, 49.244, 109.1759778, 
77.71733333, 37.446875, 101.23875, 67.78534615, 21.359, 36.54257143, 
34.13961111, 64.35253333, 80.98554545, 61.50857143, 48.983, 63.81072727, 
26.105, 46.783, 23.0605, 33.61557143, 46.31042857, 62.37061905, 
12.565, 42.31983721, 15.3982, 14.49625, 23.77425, 25.626, 74.62485714, 
170.1547778, 50.67125, 48.098, 66.83644444, 76.564875, 80.63189189, 
136.0573243, 136.3484, 86.68688889, 34.82169565, 70.00415385, 
64.67233333, 81.72766667, 57.74522034), Pass = structure(c(1L, 
2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 2L, 1L, 2L, 1L, 
2L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 
2L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 2L, 1L, 2L, 1L, 2L, 2L, 
1L, 2L, 1L, 2L), .Label = c("0", "1"), class = "factor")), .Names = c("attempt", 
"searchtime", "mean.search.flow", "Pass"), class = "data.frame", row.names = c(1L, 
2L, 3L, 4L, 5L, 6L, 7L, 8L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 
19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 28L, 29L, 30L, 31L, 32L, 
33L, 34L, 35L, 36L, 37L, 38L, 39L, 40L, 50L, 51L, 53L, 54L, 60L, 
61L, 62L, 63L, 64L, 65L, 66L, 67L, 68L, 69L, 70L, 71L, 72L))

First model with single predictor

 M2 <- glm(Pass ~ searchtime,
      data = DF3, 
      family = binomial)
 summary(M2)
 drop1(M2, test = "Chi")

Plot works fine

P1 <- predict(M2, newdata = MyData, type = "link", se = TRUE)
plot(x=MyData$searchtime, exp(P1$fit) / (1+exp(P1$fit)),
  type = "l", ylim = c(0,1),
 xlab = "search time",
 ylab = "pobability of passage")

lines(MyData$searchtime, exp(P1$fit+1.96*P1$se.fit)/
    (1 + exp(P1$fit + 1.96 * P1$se.fit)), lty = 2)
lines(MyData$searchtime, exp(P1$fit-1.96*P1$se.fit)/
    (1 + exp(P1$fit - 1.96 * P1$se.fit)), lty = 2)
points(DF3$searchtime, DF3$Search.and.pass)

Second model

M2a <- glm(Pass ~ searchtime +
        mean.search.flow+
        attempt,
      data = DF3, 
      family = binomial)
summary(M2a)
drop1(M2a, test = "Chi")

How do I plot this with "dummy" data?

I have tried along the lines of Model.matrix and expand.grid, as you would do with glmer, but fail straight away due to the two categorical variables along with factor(attempt)