would that be called a prediction interval? Maybe I'm using the wrong term)
That would not be a prediction interval. A prediction interval would incorporate uncertainty in the data generation. Its a bit useless for a binary logistic regression since we know the outcome will either be 0 or 1. A prediction interval may be more useful when you have trial data (e.g. I predict between 8 and 12 events out of the 20).
Onto calculating the confidence interval. First, start by computing $\operatorname{Var}(\operatorname{logit}(p))$. Since $\operatorname{logit}(p) = \beta_0 + \beta_1x$
$$\operatorname{Var}(\operatorname{logit}(p)) = \operatorname{Var}(\beta_0) + x^2 \operatorname{Var}(\beta_1) + 2x\operatorname{Cov}(\beta_0, \beta_1)$$
You can find these quantities from the covariance matrix estimated from the model. Take the square root of this quantity and you have the standard deviation of $\operatorname{logit}(p)$. Let's call this quantity $\sigma$ for now.
A confidence interval for $\operatorname{logit}(p)$ is $(\theta_L, \theta_U) = (\operatorname{logit}(p) - 1.96 \sigma,\operatorname{logit}(p) + 1.96 \sigma )$. Invert each endpoint to obtain a confidence interval for $p$.
Here is how one might perform this in R.
library(tidyverse)
#Simulate some data to run a regression on
x = rnorm(100)
eta = -0.8 + 0.2*x
p = plogis(eta)
y = rbinom(100, 1, p)
# Fit a model
model = glm(y~x, family = binomial())
Bigma = vcov(model)
# Formula above for standard deviation of logit p. Vectorize for ease of computation
sig =Vectorize(function(x) sqrt(Bigma[1,1] + x^2*Bigma[2,2] + 2*x*Bigma[1,2]))
# New xs to make prediction on.
new_x = seq(-2, 2, 0.01)
logit_p = predict(model, newdata=list(x=new_x))
theta_L = logit_p - 1.96*sig(new_x)
theta_U = logit_p + 1.96*sig(new_x)
# Invert the estimates
p_L = plogis(theta_L)
p = plogis(logit_p)
p_U = plogis(theta_U)
d = data.frame(x=new_x, p_L, p, p_U)
# Now compare with builtin tools
fits = predict(model, newdata = list(x=new_x), se.fit = T) %>%
bind_cols() %>%
mutate(x=new_x) %>%
mutate(p = plogis(fit),
p_L = plogis(fit - 1.96*se.fit),
p_U = plogis(fit + 1.96*se.fit))
ggplot()+
# R's computation of the confidence interval in Black
geom_line(data=fits, aes(x, p))+
geom_ribbon(data=fits, aes(x=x, ymin=p_L, ymax=p_U), alpha = 0.5)+
# Out computation of the confidence interval in Red
geom_line(data=d, aes(x, p), color = 'red')+
geom_ribbon(data=d, aes(x=x, ymin=p_L, ymax=p_U), alpha = 0.5, fill='red')
Here is an example output

The red and the black are right overtop one another, meaning the manual procedure reproduces the output of R's built in methods for calculating the standard error of the fit. You can run this code for multiple random seeds and verify they produce the same estimates (at least up to a few decimal places I imagine).