Trouble using integral in a nested f solve

I'm running an fsolve of a function which involves evaluating the integral of another function. Everything is a scalar in this problem. Every line prior to the integral line runs fine and spits out a scalar answer. However when I try to evaluate the integral it says "Subscripted assignment dimension mismatch". Any help would be appreciated! Here's the code: I'm running it as a Monte Carlo but we can just set all parameter and variable values: gamma=-1; beta=1; c=7;

    x1= 1; 
    x2= 1; 
    xi1=0.2;  
    xi2=0.21;  
    mc1=7;
    mc2=7;
    xxi=zeros(1,2); 
    xxi0 = [0,0]; 
    pd0 = [0,0]; 
pm1=8;
pm2=8;
    pm0=0;

The main function:

function F = func(xxi, gamma, beta, mc1, mc2, pm1, pm2, x1, x2, pd1, pd2, xi1, xi2)

xxi1 = xxi(1); xxi2 = xxi(2);

fun3 = @(a, b)funpd21(gamma, beta, mc1, mc2, x1, x2, a, b); fun4 = @(c, d)funpd22(gamma, beta, mc1, mc2, x1, x2, c, d);

fun1 = @(xxi1,e)(fun3(xxi1,e)-mc1).*exp(x1*beta+fun3(xxi1,e).*gamma+xxi1)./(1+exp(x1.*beta+fun3(xxi1,e).*gamma+xxi1)+exp(x2.*beta+fun4(xxi1,e).*gamma+e)).*normpdf(e,0.2,1); fun2 = @(f,xxi2)(fun4(f,xxi2)-mc2).*exp(x2*beta+fun4(f,xxi2).*gamma+xxi2)./(1+exp(x1.*beta+fun3(f,xxi2).*gamma+f)+exp(x2.*beta+fun4(f,xxi2).*gamma+xxi2)).*normpdf(f,0.2,1);

F(1) = normcdf(xxi(2),0.2,1).*(pm1-mc1).*exp(x1*beta+pm1.*gamma+xxi1)./(1+exp(x1.*beta+pm1.*gamma+xxi1))+(1-normcdf(xxi(2),0.2,1))*integral(@(g)fun1(xxi1,g),xxi2,5)/(1-normcdf(xxi2,0.2,1)); F(2) = normcdf(xxi(1),0.2,1).*(pm2-mc2).*exp(x2*beta+pm2.*gamma+xxi2)./(1+exp(x2.*beta+pm2.*gamma+xxi2))+(1-normcdf(xxi(1),0.2,1))*integral(@(h)fun2(h,xxi2),xxi1,5)/(1-normcdf(xxi1,0.2,1)); end

And the functions earlier in the nest: function pd91 = funpd21(gamma, beta, mc1, mc2, x1,x2,xi1,xi2)

options = optimset('Display','off');

f21 = @(pd)funpd(pd, gamma, beta, mc1, mc2, x1,x2,xi1,xi2); pd9 = fsolve(f21, [0,0], options);

pd91 = pd9(1);

end

function pd92 = funpd22(gamma, beta, mc1, mc2, x1,x2,xi1,xi2)

options = optimset('Display','off');

f20 = @(pd)funpd(pd, gamma, beta, mc1, mc2, x1,x2,xi1,xi2); pd9 = fsolve(f20, [0,0], options);

pd92 = pd9(2);

end

function F = funpd(pd, gamma, beta, mc1, mc2, x1, x2, xi1, xi2) pd1 = pd(1); pd2 = pd(2); F(1) = (pd1-mc1).*gamma+1./(1-exp(x1.*beta+pd1.*gamma+xi1)./(1+exp(x1.*beta+pd1.*gamma+xi1)+exp(x2.*beta+pd2.*gamma+xi2))); F(2) = (pd2-mc2).*gamma+1./(1-exp(x2.*beta+pd2.*gamma+xi2)./(1+exp(x1.*beta+pd1.*gamma+xi1)+exp(x2.*beta+pd2.*gamma+xi2))); end