2d double integral with trapz

Hi: I am using Comsol to connect to matlab and run the following code. For some reason, the function trapz crashes. It was working before. Can any one help me with this?. what is wrong with this function ?. Was there any change in the use of this function in recent versions of matlab? please help. Naema

%Set parameters to compute fieber field.
lambda0 = 830e-9;
x0 = 0E-6; y0 = 0E-6;       %x and y fiber offset
%%a - core radius;
%for SMF - 28: a = 4.1e-6; NA = .12 @ 1550nm;
a = 2e-6; NA = .12; V = 2*pi*a/lambda0*NA;
r_f = (5.3/2)*1E-6; %a*(.65 + 1.619/V^(3/2) + 2.879/V^6);

%Extract SPP field from Femlab.
fem.xmesh=meshextend(fem);
gridNumIntervals = 800;
Ymax = 5e-6;
Xmax = 5e-6;
[x, y] = meshgrid(-Xmax:2*Xmax/gridNumIntervals:Xmax, -Ymax:2*Ymax/gridNumIntervals:Ymax);
p=[x(:)'; y(:)'];
Espp=postinterp(fem,'Ey',p);
Espp=reshape(Espp, size(x));

Efib = exp(-((x-x0).^2 + (y-y0).^2)./r_f^2);

if 0
    %Some times Femlab does not want to plot, so better save the fields and
    %plot them in matlab.

      figure(1);
      subplot(211);
      hold on;
      surf(x*1e6, y*1e6, abs(Espp));
      colorbar;
      shading flat;

      axis image;
      xlabel('x(\mum)', 'fontsize', 12);
      ylabel('y(\mum)', 'fontsize', 12);
      % colormap(gray) %good colours 'bone' 'cool' 'winter'

      subplot(212);
      hold on;
      surf(x*1e6, y*1e6, abs(Efib));
      colorbar;
      shading flat;

      axis image;
      xlabel('x(\mum)', 'fontsize', 12);
      ylabel('y(\mum)', 'fontsize', 12);
      %colormap(gray) %good colours 'bone' 'cool' 'winter'

end;

%Save the SPP and fiber fields and the coordinates.
% save('C:\Working folder\Theoretical\spp coupling eff\CEFF.mat','x','y','Espp','Efib');

save('C:\kholoud desk_2\Haneen\Israel prob\CEFF.mat','x','y','Espp','Efib');

%Compute coupling efficeincy.
fn=Efib.*conj(Espp);
fd=sqrt(trapz2D(x,y,Espp.*conj(Espp))*trapz2D(x,y,Efib.*conj(Efib)));
C = trapz2D(x,y,fn./fd);
absC = abs(C);
Ceff = absC;
CdB = -10*log10(Ceff);

beta0=2*pi/lambda0;
beta = sqrt(-fem.sol.lambda); %from the fem.sol.lambda
neff = beta/beta0; %beta*lambdaSweep(fInt)/2/pi;
MPA=20*log10(exp(1))*imag(beta/100)

%Display results.
disp('--------------------------------------------');
disp('');
disp(['Grid size in x: ', num2str(2*Xmax*1E6/gridNumIntervals),' um']);
disp(['Grid size in y: ', num2str(2*Ymax*1E6/gridNumIntervals),' um']);
disp('');
disp(['Fiber offset in x:', num2str(x0*1E6),' um']);
disp(['Fiber offset in y:', num2str(y0*1E6),' um']);
disp(['Complex coupling factor: ', num2str(C)]);
disp(['Coupling efficiency: ', num2str(Ceff)]);
disp(['Coupling loss per facet: ', num2str(CdB), ' dB']);

##########################################

here is the function associated with the code: % % A double integration of the tabulated data. % f(x,y) must be genarated in terms of meshgrid % Integral computes by approximation of the integral of f(x,y) via the trapezoidal method. % Firstly, an integration loop within inner integral over X is % performed obtaining the function of area A(y). Secondly, outer loop % integration is performed for A(y) over y limits % % Example: % Q = dblquad(inline('y*sin(x)+x*cos(y)'), pi, 2*pi, 0, pi); % Q = -9.8696 % [X,Y] = meshgrid(pi:pi/100:2*pi, 0:pi/100:pi ); % Z=Y.*sin(X)+X.*cos(Y); % int_2D_tabulated(X,Y, Z ) % ans = -9.8688

function V = trapz2D(X, Y, Z )

V = trapz( Y(:,1), (trapz(X(1,:),Z, 2)) ); return