index exceeds matrix dimensions

Index exceeds matrix dimensions.

Error in
HMWK2_P3>@(theta1,phi)((1./2.*pi()).*((pi()-phi).*cos(phi)+sin(phi)).*tan(theta0(i)).*tan(theta1)-(1./pi()).*(tan(theta0(i))+tan(theta1)+sqrt(tan(theta1).^2+tan(theta0(i)).^2-2.*tan(theta0(i)).*tan(theta1).*cos(phi)))).*cos(theta1(i)).*sin(theta1)

Error in integral2Calc>@(y)fun(xi*ones(size(y)),y) (line 18)
    @(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions), ...

Error in integralCalc/iterateScalarValued (line 323)
                fx = FUN(t).*w;

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in
integral2Calc>@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions)

Error in
integral2Calc>@(x)arrayfun(@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions),x,ymin(x),ymax(x))
(line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( ...

Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral2Calc>integral2i (line 20)
[q,errbnd] = integralCalc(innerintegral,xmin,xmax,opstruct.integralOptions);

Error in integral2Calc (line 7)
    [q,errbnd] = integral2i(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 106)
    Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

Error in HMWK2_P3 (line 20)
I1cs(i) = integral2
(fun1,0,pi/2,0,pi,'Method','iterated','AbsTol',0,'RelTol',1e-10);

The code is:

    %%Homework 2 Problem 3
    clear;
    clc;
    theta0=zeros(1,76); %Preallocation
    I1=zeros(1,76);
    I1cs=zeros(1,76);
    I2=zeros(1,76);
    I2cs=zeros(1,76);

    for i = 1:76
        theta0(i)=deg2rad(i-1);

    %z=(@(theta1,phi) (acos(cos(theta1).*cos(theta0(i))+sin(theta1).*sin(theta0(i)).*cos(phi)));

    % Given functions f1 and f2
    fun1= (@(theta1,phi) ((1./2.*pi()).*((pi()-phi).*cos(phi)+sin(phi)).*tan(theta0(i)).*tan(theta1)-(1./pi()).*(tan(theta0(i))+tan(theta1)+sqrt(tan(theta1).^2 + tan(theta0(i)).^2 - 2.* tan(theta0(i)).*tan(theta1).*cos(phi)))).*cos(theta1(i)).*sin(theta1));
    fun2= (@(theta1,phi) ((4./(3.*pi())).*(1./(cos(theta1)+cos(theta0(i)))).*(((pi()./2)-(acos(cos(theta1).*cos(theta0(i))+sin(theta1).*sin(theta0(i)).*cos(phi)))).*cos((acos(cos(theta1).*cos(theta0(i))+sin(theta1).*sin(theta0(i)).*cos(phi))))+sin((acos(cos(theta1).*cos(theta0(i))+sin(theta1).*sin(theta0(i)).*cos(phi)))))-(1./3)).*cos(theta1(i)).*sin(theta1));

    I1cs(i) = integral2  (fun1,0,pi/2,0,pi,'Method','iterated','AbsTol',0,'RelTol',1e-10);
    I2cs(i)= integral2 (fun2,0,pi/2,0,pi,'Method','iterated','AbsTol',0,'RelTol',1e-10);
    I1(i)= -0.9946 - 0.0281 * tan(theta0(i)) - 0.0916 * tan(theta0(i))^2 +0.0108 * tan(theta0(i))^3;
    I2(i)= -0.0137 + 0.0370 * tan(theta0(i)) + 0.0310 * tan(theta0(i))^2 - 0.0059 * tan(theta0(i))^3;

    end

I can't figure out where the disagreement is.