Solving a system of two integral equations numerically

I am trying to solve a system of two equations for given that and letting , vary from 0.1 to 1 in steps of 0.1.

I think I need to do a double for loop that numerically approximates the integral and then solves for appropriate values.

I did make an attempt but it is not going so well.

count= 1; %counter
svals1 = zeros(10,20); %an array
svals2 = zeros(10,20); %an array
g = 0.65; %gamma
d1= 1; %delta1
b1 = 1.1; %beta1
a11= 1; %alpha_11
a12 = 1; %alpha_12
a21 = 1; %alpha_21
a22 = 1; %alpha_22
syms x;
syms S1;
syms S2;

for d2= 0:0.1:count
    for b2 = 0:0.1:count
      svals1(d2,b2,x,S1,S2) = vpasolve( (b1.*g)/d1.*integral( exp(- g.*x).* (1-exp(-d1.* x)).* exp(-(a11.*(b1.*S1)/d1 (1-exp(-d1.* x) ) +a12.* (b2.*S2)/d2 (1-exp(-d2.* x)))),0,Inf)-1,S1)
      svals2(d2,b2,x,S1,S2) = vpasolve( (b2.*g)/d2.*integral( exp(- g.* x).* (1-exp(-d2.* x)) .*exp(-(a22.*(b2.*S2)/d2 (1-exp(-d2.* x) ) +a21.* (b1.*S1)/d1 (1-exp(-d1.* x)))),0,Inf)-1,S2)
      disp(d2,b2,svals1,svals2)
    end
end