How to draw posterior samples from the model with multiple outcomes?

Hi, my model is as follows :

model =
bf(Y1 ~ IV1+…IVn+(1|rand|individual)+factor(time), family=bernoulli(logit)) +
bf(Y2 ~ IV1+…IVn+(1|rand|individual)+factor(time), family=bernoulli(logit)) +
bf(DV1 | subset(Y1) ~ IV1+…IVn+(1|rand|individual)+factor(time), family=gaussian()) +
bf(DV2 | subset(Y1) ~ IV1+…IVn+(1|rand|individual)+factor(time), family=gaussian())
Prior_int = c(
set_prior(“normal(0,5)”, resp = "Y1, class = “Intercept”),
set_prior(“normal(0,5)”, resp = “Y2”, class = “Intercept”),
set_prior(“normal(0,5)”, resp = “DV1”, class = “Intercept”),
set_prior(“normal(0,5)”, resp = “DV2”, class = “Intercept”) )
fit = brm (model + set_rescor(FALSE), data=DATASET,
warmup = 2000, iter = 4000, seed = 2021,
chains = 4, control = list(adapt_delta =.80),
cores = 4, prior = c(Prior_int), inits = 0)

What I want to do is to obtain posterior samples from different values of IV1 (e.g., IV1=0, IV1=50, IV1=100) with each of the outcomes (e.g., Y1,Y2,…) and put them in a plot to compare the fitted values.
I found that I could get population-level (fixed effects) posterior draws using “fixef(fit, summary=F)”.

And I want to know

  1. How could I get 95% credible interval when calculating the fitted values using fixef for all of the draws (as I know, it is not possible to use the ‘prob’ option when ‘summary=F’ is designated).
  2. I saw a posting that recommends to use the function ‘fitted’ , but when I use the function error comes out (see below)
    > fitted(fit_main, newdata=data.frame(IV1=xx), scale = “response”, summary=FALSE)
    Error: Argument ‘resp’ must be a single variable name for models using addition argument ‘draw_ids’.
  3. how to draw a plot for comparing the fitted value of two outcomes (e.g., comparing the fitted value of 8,000 draws for Y1 and Y2 when IV1=0 in the same plot)

Thank you in advance for your help.