MATLAB: 2-D plot with z-axis given in color

My friends and I have been struggling to generate a 2-D plot in MATLAB with eta_1 and eta_2 both varying in 0:0.01:1 and the z-axis given by color.

We have a system of 8 differential equations, with HIVinf representing the total new HIV infections in a population over 1 year (HIVinf is obtained by integrating a function of eta_1, eta_2).

We are looping through eta_1 and \eta_2 (two 'for' loops) with the ode45 solver within the 'for' loops.

Based on our prior numerical results, we should be getting much color variation in the 2D-plot. There should be patterns of darkness (high concentration of HIVinfections) along the edges of the plot, and lightness along the diagonals (low concentrations).

However, the following snippet does not produce what we want (I have attached the figure below).

    [X,Y] = meshgrid(eta_11,eta_22);
    figure;
    pcolor(X,Y,AA); 
    shading interp;

[1]: http://i.stack.imgur.com/XNM7b.jpg

I have attached the code below, as concisely as possible. The function ydot works fine (it is required to run ode45).

*We would greatly appreciate if you could help us fix the snippet.*

----------

    function All()

    global Lambda mu mu_A mu_T beta tau eta_1 eta_2 lambda_T rho_1 rho_2 gamma
    alpha = 20;
    TIME = 365;
    eta_11 = zeros(1,alpha);
    eta_22 = zeros(1,alpha); 
    AA = zeros(1,alpha);
    BB = zeros(1,alpha);
    CC = zeros(1,alpha);

    for n = 1:1:alpha
        for m = 1:1:alpha 
        Lambda = 531062; 
        mu = 1/25550;
        mu_A = 1/1460;
        mu_T = 1/1825;
        beta = 187/365000; 
        tau = 4/365;  
        lambda_T = 1/10;
        rho_1 = 1/180;
        rho_2 = 1/90;
        gamma = 1/1000;

        eta_1 = (n-1)./(alpha-1);
        eta_11(m) = (m-1)./(alpha-1);

        eta_2 = (m-1)./(alpha-1);
        eta_22(m) = (m-1)./(alpha-1);

        y0 = [191564208, 131533276, 2405629, 1805024, 1000000, 1000000, 500000, 500000];

        [t,y] = ode45('SimplifiedEqns',[0:1:TIME],y0);

        N = y(:,1)+y(:,2)+y(:,3)+y(:,4)+y(:,5)+y(;,6)+y(:,7)+y(:,8);

        HIVinf1=[0:1:TIME];
        HIVinf2=[beta.*(S+T).*(C1+C2)./N];
        HIVinf=trapz(HIVinf1,HIVinf2);
        AA(n,m) = HIVinf;
        end    
    end

    [X,Y] = meshgrid(eta_11,eta_22);
    figure;
    pcolor(X,Y,AA); 
    shading interp;

    function ydot = SimplifiedEqns(t,y)

    global Lambda mu mu_A mu_T beta tau eta_1 eta_2 lambda_T rho_1 rho_2 gamma

    S = y(1);
    T = y(2);
    H = y(3); 
    C = y(4);
    C1 = y(5); 
    C2 = y(6);
    CM1 = y(7); 
    CM2 = y(8);

    N = S + T + H + C + C1 + C2 + CM1 + CM2;
    ydot = zeros(8,1);

    ydot(1)=Lambda-mu.*S-beta.*(H+C+C1+C2).*(S./N)-tau.*(T+C).*(S./N);
    ydot(2)=tau.*(T+C).*(S./N)-beta.*(H+C+C1+C2).*(T./N)-(mu+mu_T).*T;
    ydot(3)=beta.*(H+C+C1+C2).*(S./N)-tau.*(T+C).*(H./N)-(mu+mu_A).*H;
    ydot(4)=beta.*(H+C+C1+C2).*(T./N)+tau.*(T+C).*(H./N)-   (mu+mu_A+mu_T+lambda_T).*C;
    ydot(5)=lambda_T.*C-(mu+mu_A+rho_1+eta_1).*C1;
    ydot(6)=rho_1.*C1-(mu+mu_A+rho_2+eta_2).*C2;
    ydot(7)=eta_1.*C1-(mu+rho_1+gamma).*CM1;
    ydot(8)=eta_2.*C2-(mu+rho_2+gamma.*(rho_1)./(rho_1+rho_2)).*CM2+(rho_1).*CM1;
    end
    end