Double Integral with Varying Limits

syms x y
t=(30:1:90);
nu=0.3;
ms=1+nu;
E=2.07*10^11;
G=8.28*10^10;
R=0.0185;
Ep=E/(1-nu^2);
k=(6*(1+nu))/(7+6*nu);
A=pi*R^2;
L=0.254;
sigma=1;
I=(1/4)*pi*R^4;
mu=R*x;
mup=2*R*(1-sigma+x);
beta=R*y;
for i=1:length(t)
F1=(sqrt(((2*mup)/(pi*mu))*tan((pi*mu)/(2*mup))))*((0.752+2.02*(mu/mup)+0.37*(1-sin((pi*mu)/(2*mup)))^3)/(cos((pi*mu)/(2*mup))));
Flll=sqrt(((2*mup)/(pi*mu))*tan((pi*mu)/(2*mup)));
x_s=(((0.5*L)/R)+y*cosd(t(i)));
  f11=matlabFunction(32*(x_s.^2)*(y^2)*x*(((sind(t(i))).^6)*F1^2+ms*((sind(t(i))).^4).*((cosd(t(i))).^2)*Flll^2)+...
      16*x_s*k*y*x*((sind(t(i)).^3).*sind(2*t(i))*F1^2+ms*(sind(t(i)).^2).*cosd(t(i)).*cosd(2*t(i))*Flll^2)+...
      2*k^2*x*(F1^2*(sind(2*t(i)).^2)+ms*(cosd(2*t(i))).^2*Flll^2));
  xmin=0;
  xmax=(sqrt(1-y^2*(sind(t(i))).^2))-1+sigma;
  ymax=(sqrt(1-(1-sigma)^2))/(sind(t(i)));
  ymin=-ymax;
  Q(i)=(L)/(k*G*A)+(L^3)/(3*E*I)+1/(Ep*pi*R)*int(int(f11,x,xmin,xmax),y,ymin,ymax); % double integral
end
plot(t,Q,'g','LineWidth',3);grid on;
title('Flexibility Coefficient');
xlabel('theta');
ylabel('f11');
The above code that I write it, but when I run for sigma=0.5 don't response.
How to plot the solution for sigma=0:0.01:1 and t=20:1:180 ?
best regards