Solving a system of integral equations in R using multiroot

I have a system of integrals that have free variables S1 and S2. I would like to know what values of S1, S2 make the resulting (integrals-1) equal to 0. (I will eventually change the parameter values so that they are not the same equations.)

Multiroot uses the library: library(rootSolve).

a11 <- 1 #alpha_{11}
a12 <- 1 #alpha_{12}
a21 <- 1 #alpha_{21}
a22 <- 1 #alpha_{22}
b1 <- 2  #beta1
b2 <- 2 #beta2
d1 <- 1 #delta1
d2 <- 1 #delta2
g <- 0.5 #gamma


integrand1 <- function(x,S) {b1*g/d1*exp(-g*x)*(1-exp(-d1* x))*exp(-a11*b1*S[1]/d1*(1-exp(-d1*x))-a12*b2*S[2]/d2*(1-exp(-d2*x)))}
integrand2 <- function(x,S) {b2*g/d2*exp(-g*x)*(1-exp(-d2* x))*exp(-a22*b2*S[2]/d2*(1-exp(-d2*x))-a21*b1*S[1]/d1*(1-exp(-d1*x)))}

#defining equation we would like to solve
intfun1<- function(S) {integrate(function(x)integrand1(x, S),lower=0,upper=Inf)[[1]]-1}
intfun2<- function(S) {integrate(function(x)integrand2(x, S),lower=0,upper=Inf)[[1]]-1}

#putting both equations into one term
model <- function(S) c(F1 = intfun1,F2 = intfun2)

#Solving for roots
(ss <-multiroot(f=model, start=c(0,0)))

Gives the error:

Error in stode(y, times, func, parms = parms, ...) : REAL() can only be applied to a 'numeric', not a 'list'

For reference this simpler system works fine:

#Defining the functions
model <- function(x) c(F1 = x[1]+ 4*x[2] -8,F2 = x[1]-4*x[2])

#Solving for the roots
(ss <- multiroot(f = model, start = c(0,0)))

If I am mistaken in using multiroot for solving the roots of a system of equations and you think there is a better function I am open to suggestions. Thank you!