How to speed up the program?

(x,y) is a function generated as the following codes. I need to calculate the integration that "integral2(@(X,Y) arrayfun(@(x,y) Tuu(x,y),X,Y),0.01,pi/2,0,pi/4)". However, the program is too slow to use. How to speed up this program in details? Many thanks! The codes of Tuu(x,y) are as following:

function U=Tuu(x,y)
  syms kz d
  m = 2;
  dd=2.106*(m+1);
  vh = 4;
  mu = 11;
  delta = 8;
  HBAR = 1.05457266e-34;
  ME = 9.1093897e-31;
  ELEC = 1.60217733e-19;
  Kh = 2.106;
  vKh = [0,0,0;Kh,0,0;-Kh,0,0;0,Kh,0;0,-Kh,0];
  kc = sqrt(2.*ME.*ELEC/HBAR^2).*1e-10;
  ku = kc.*sqrt(mu+delta);
  kd = kc.*sqrt(mu-delta);
  a3 = [pi/Kh,pi/Kh,sqrt(2).*pi/Kh];
  kuu =[-ku.*sin(x).*cos(y), -ku.*sin(x).*sin(y), kz];
  n=0:m;
  for p=1:5;
      for q=1:5;
      tuu(p,q)= (sum((kuu + vKh(p,:)).^2)-ku^2).*(p==q)+ kc^2*vh*sum(exp(i.*n.*sum((vKh(q,:)-vKh(p,:)).*a3)))/(m+1).*(p~=q);
      end
  end
  dtuu=det(tuu);
  kz0=vpasolve(dtuu,kz);
  kzz=kz0(real(vpa(kz0))>=0&imag(vpa(kz0))>=0);
  tuu1=subs(tuu,kz,kzz(1));
  tuu2=subs(tuu,kz,kzz(2));
  tuu3=subs(tuu,kz,kzz(3));
  tuu4=subs(tuu,kz,kzz(4));
  tuu5=subs(tuu,kz,kzz(5));
  tuu11=double(tuu1);
  tuu22=double(tuu2);
  tuu33=double(tuu3);
  tuu44=double(tuu4);
  tuu55=double(tuu5);
  nuu1=null(tuu11);
  nuu2=null(tuu22);
  nuu3=null(tuu33);
  nuu4=null(tuu44);
  nuu5=null(tuu55);
  piuu=[nuu1,nuu2,nuu3,nuu4,nuu5];
  pei=[1;0;0;0;0];
  A=piuu\pei;
  psiuu1=A(1).*nuu1(1)*exp(i*kzz(1)*d)+A(2)*nuu2(1)*exp(i*kzz(2)*d)+A(3)*nuu3(1)*exp(i*kzz(3)*d)+A(4)*nuu4(1)*exp(i*kzz(4)*d)+A(5)*nuu5(1)*exp(i*kzz(5)*d);
  psiuu2=A(1).*nuu1(2)*exp(i*kzz(1)*d)+A(2)*nuu2(2)*exp(i*kzz(2)*d)+A(3)*nuu3(2)*exp(i*kzz(3)*d)+A(4)*nuu4(2)*exp(i*kzz(4)*d)+A(5)*nuu5(2)*exp(i*kzz(5)*d);
  psiuu3=A(1).*nuu1(3)*exp(i*kzz(1)*d)+A(2)*nuu2(3)*exp(i*kzz(2)*d)+A(3)*nuu3(3)*exp(i*kzz(3)*d)+A(4)*nuu4(3)*exp(i*kzz(4)*d)+A(5)*nuu5(3)*exp(i*kzz(5)*d);
  psiuu4=A(1).*nuu1(4)*exp(i*kzz(1)*d)+A(2)*nuu2(4)*exp(i*kzz(2)*d)+A(3)*nuu3(4)*exp(i*kzz(3)*d)+A(4)*nuu4(4)*exp(i*kzz(4)*d)+A(5)*nuu5(4)*exp(i*kzz(5)*d);
  psiuu5=A(1).*nuu1(5)*exp(i*kzz(1)*d)+A(2)*nuu2(5)*exp(i*kzz(2)*d)+A(3)*nuu3(5)*exp(i*kzz(3)*d)+A(4)*nuu4(5)*exp(i*kzz(4)*d)+A(5)*nuu5(5)*exp(i*kzz(5)*d);
  Tuux=ku.*sin(x).*(imag(conj(psiuu1)*diff(psiuu1,d)+conj(psiuu2)*diff(psiuu2,d)+conj(psiuu3)*diff(psiuu3,d)+conj(psiuu4)*diff(psiuu4,d)+conj(psiuu5)*diff(psiuu5,d)));
  U=double(subs(Tuux,d,dd));
  end

[Stephen23]

" How to speed up this program in details?"

Probably by performing numeric calculations rather than symbolic ones.

[Henan Fang]

@Stephen Cobeldick Can you tell me, in detail, how to perform numeric calculations instead of symbolic ones in this program?

[Stephen23]

"Can you tell me, in detail, how to perform numeric calculations instead of symbolic ones in this program?"

Not really, as that function contains no code comments, no help, no links, no references, no examples, no demo, no explanation what it is doing. But it is possible that you could replace symbolic operations with faster numeric operations, if you want a numeric result at the end.

If you ask the author of that function they should be able to describe what the function does, and then you can start to investigate numeric alternatives.