title: nested trapz double integration
body: <p>I'd like knowing if there is any way to bypass the call for <code>quad2d</code> with a nested <code>trapz</code> loop. I'll discuss my problem with some more detail: currently, I perform the calculation of a double integral this way:</p>
  
  <pre><code>clc, clear all, close all
  load E_integral.mat
  
  c = 1.476;
  gamma = 3.0;
  
  beta_int = (c*gamma)./(k_int.*sqrt(E_integral));
  
  figure, loglog(k_int,beta_int,'r','LineWidth',2.0), grid on;
  
  k1 = (.01:.1:100);
  k2 = .01:.1:100;
  k3 = -100:.1:100;
  
  int_k3 = zeros(size(k2));
  int_k3k2 = zeros(size(k1));
  
  tic
  for ii = 1:numel(k1)
      phi11 = @(k2,k3) PHI11(k1(ii),k2,k3,k_int,beta_int);
      int_k3(ii) = 2*quad2d(phi11,-100,100,-100,100);
  end
  toc
  </code></pre>
  
  <p>where <code>PHI11</code> is defined as</p>
  
  <pre><code>function phi11 = PHI11(k1,k2,k3,k_int,beta_int)
  k = sqrt(k1.^2 + k2.^2 + k3.^2);
  ksq = k.^2;
  k1sq = k1.^2;
  fourpi = 4.*pi;
  beta = exp(interp1(log(k_int),log(beta_int),log(k),'linear'));
  k30 = k3 + beta.*k1;
  k0 = sqrt(k1.^2 + k2.^2 + k30.^2);
  k0sq = k0.^2;
  k04sq = k0.^4;
  Ek0 = (1.453.*k04sq)./((1 + k0sq).^(17/6));
  
  C1 = (beta.*k1sq.*(k0sq - 2.*(k30.^2) + beta.*k1.*k30))./(ksq.*(k1.^2 + k2.^2));
  C2 = ((k2.*k0sq)./((k1.^2 + k2.^2).^(3/2))).*atan2((beta.*k1.*sqrt(k1.^2 + k2.^2)),(k0sq - k30.*k1.*beta));
  xhsi1 = C1 - (k2./k1).*C2;
  xhsi1_sq = xhsi1.^2;
  phi11 = (Ek0./(fourpi.*k04sq)).*(k0sq - k1sq - 2.*k1.*k30.*xhsi1 + (k1.^2 + k2.^2).*xhsi1_sq);
  end
  </code></pre>
  
  <p>and <code>E_integral.mat</code> can be obtained this way:</p>
  
  <pre><code>clc,clear all,close all
  
  k_int = .001:.01:1000;
  
  Ek = (1.453.*k_int.^4)./((1 + k_int.^2).^(17/6));
  
  
  E = @(k_int) (1.453.*k_int.^4)./((1 + k_int.^2).^(17/6));
  
  E_integral = zeros(size(k_int));
  
  for ii = 1:numel(k_int)    
      E_integral(ii) = integral(E,k_int(ii),Inf);    
  end
  
  save('E_integral','k_int','E_integral')
  </code></pre>
  
  <p>Now the question is: would it be possible to overlook <code>quad2d</code> and the <code>handle function</code> in favor of a more practicle approach, by using a nested <code>trapz</code> function?</p>
  
  <p>So far, I've tried the following piece of code, which has not yielded to the expected results:</p>
  
  <pre><code>int_k33 = zeros(size(k2));
  S_11 = zeros(size(k1));
  tic
  for ii = 1:1
      for jj = 1:numel(k2)
          int_k33(jj) = trapz(k3,PHI11(k1(ii),k2(jj),k3,k_int,beta_int));        
      end
      S_11(ii) = 4*trapz(k2,int_k33);
  end
  toc
  </code></pre>