Stepwisefit changing yy variable...

I have written a code for a stepwise regression. Usually there would be one yy variable and multiple xx variables. Part of the difference between my xx variables is a certain lag (0-6).

The problem is that where a lag exists in an xx variable, it must be matched in the yy variable, so the yy variable needs to change according to the xx variable and vice versa. See the code below (Not including loading data). Thank you in advance for any help

        for j = 1:1505
          Nc0=NDVI(1:end,j);
          A1c0=Approx1 (1:end-0,j);
          A2c0=Approx2 (1:end-0,j);  
          A3c0=Approx3 (1:end-0,j);
          A4c0=Approx4 (1:end-0,j);
          D1c0=Det1 (1:end-0,j);
          D2c0=Det2 (1:end-0,j);
          D3c0=Det3 (1:end-0,j);
          D4c0=Det4 (1:end-0,j);
          Nc1=NDVI(1+1:end,j);
          A1c1=Approx1 (1:end-1,j);
          A2c1=Approx2 (1:end-1,j);  
          A3c1=Approx3 (1:end-1,j);
          A4c1=Approx4 (1:end-1,j);
          D1c1=Det1 (1:end-1,j);
          D2c1=Det2 (1:end-1,j);
          D3c1=Det3 (1:end-1,j);
          D4c1=Det4 (1:end-1,j);
          Nc2=NDVI(1+2:end,j);
          A1c2=Approx1 (1:end-2,j);
          A2c2=Approx2 (1:end-2,j);  
          A3c2=Approx3 (1:end-2,j);
          A4c2=Approx4 (1:end-2,j);
          D1c2=Det1 (1:end-2,j);
          D2c2=Det2 (1:end-2,j);
          D3c2=Det3 (1:end-2,j);
          D4c2=Det4 (1:end-2,j);
          Nc3=NDVI(1+3:end,j);
          A1c3=Approx1 (1:end-3,j);
          A2c3=Approx2 (1:end-3,j);  
          A3c3=Approx3 (1:end-3,j);
          A4c3=Approx4 (1:end-3,j);
          D1c3=Det1 (1:end-3,j);
          D2c3=Det2 (1:end-3,j);
          D3c3=Det3 (1:end-3,j);
          D4c3=Det4 (1:end-3,j);
          Nc4=NDVI(1+4:end,j);
          A1c4=Approx1 (1:end-4,j);
          A2c4=Approx2 (1:end-4,j);  
          A3c4=Approx3 (1:end-4,j);
          A4c4=Approx4 (1:end-4,j);
          D1c4=Det1 (1:end-4,j);
          D2c4=Det2 (1:end-4,j);
          D3c4=Det3 (1:end-4,j);
          D4c4=Det4 (1:end-4,j);
          Nc5=NDVI(1+5:end,j);
          A1c5=Approx1 (1:end-5,j);
          A2c5=Approx2 (1:end-5,j);  
          A3c5=Approx3 (1:end-5,j);
          A4c5=Approx4 (1:end-5,j);
          D1c5=Det1 (1:end-5,j);
          D2c5=Det2 (1:end-5,j);
          D3c5=Det3 (1:end-5,j);
          D4c5=Det4 (1:end-5,j);
          Nc6=NDVI(1+0:end,j);
          A1c6=Approx1 (1:end-6,j);
          A2c6=Approx2 (1:end-6,j);  
          A3c6=Approx3 (1:end-6,j);
          A4c6=Approx4 (1:end-6,j);
          D1c6=Det1 (1:end-6,j);
          D2c6=Det2 (1:end-6,j);
          D3c6=Det3 (1:end-6,j);
          D4c6=Det4 (1:end-6,j);
          yy = Nc0; %THIS IS WHERE THE PROBLEM OCCURS, I NEED yy TO INCLUDE Nc0, Nc1, Nc2 ETC BUT ONLY TO %WORK WITH THE CORRESPONDING xx
          xx=[A1c0, A2c0, A3c0, A4c0, D1c0, D2c0, D3c0, D4c0,A1c1, A2c1, A3c1, A4c1, D1c1, D2c1, D3c1, D4c1,...
              A1c2, A2c2, A3c2, A4c2, D1c2, D2c2, D3c2, D4c2,A1c3, A2c3, A3c3, A4c3, D1c3, D2c3, D3c3, D4c3,...
              A1c4, A2c4, A3c4, A4c4, D1c4, D2c4, D3c4, D4c4,A1c5, A2c5, A3c5, A4c5, D1c5, D2c5, D3c5, D4c5,...
              A1c6, A2c6, A3c6, A4c6, D1c6, D2c6, D3c6, D4c6];
          %Begin Stepwise Regression
          if isnan(Nc0)
               continue
           else
                  [B,SE,PVAL,INMODEL,STATS,NEXTSTEP,HISTORY]= ...
                                           stepwisefit(xx,yy,'penter',.05);
            A1nolag(j)=INMODEL(1);
            A2nolag(j)=INMODEL(2);
            A3nolag(j)=INMODEL(3);
            A4nolag(j)=INMODEL(4);
            D1nolag(j)=INMODEL(5);
            D2nolag(j)=INMODEL(6);
            D3nolag(j)=INMODEL(7);
            D4nolag(j)=INMODEL(8);
            A11lag(j)=INMODEL(9);
            A21lag(j)=INMODEL(10);
            A31lag(j)=INMODEL(11);
            A41lag(j)=INMODEL(12);
            D11lag(j)=INMODEL(13);
            D21lag(j)=INMODEL(14);
            D31lag(j)=INMODEL(15);
            D41lag(j)=INMODEL(16);
            A12lag(j)=INMODEL(17);
            A22lag(j)=INMODEL(18);
            A32lag(j)=INMODEL(19);
            A42lag(j)=INMODEL(20);
            D12lag(j)=INMODEL(21);
            D22lag(j)=INMODEL(22);
            D32lag(j)=INMODEL(23);
            D42lag(j)=INMODEL(24);
            A13lag(j)=INMODEL(25);
            A23lag(j)=INMODEL(26);
            A33lag(j)=INMODEL(27);
            A43lag(j)=INMODEL(28);
            D13lag(j)=INMODEL(29);
            D23lag(j)=INMODEL(30);
            D33lag(j)=INMODEL(31);
            D43lag(j)=INMODEL(32);
            A14lag(j)=INMODEL(33);
            A24lag(j)=INMODEL(34);
            A34lag(j)=INMODEL(35);
            A44lag(j)=INMODEL(36);
            D14lag(j)=INMODEL(37);
            D24lag(j)=INMODEL(38);
            D34lag(j)=INMODEL(39);
            D44lag(j)=INMODEL(40);
            A15lag(j)=INMODEL(41);
            A25lag(j)=INMODEL(42);
            A35lag(j)=INMODEL(43);
            A45lag(j)=INMODEL(44);
            D15lag(j)=INMODEL(45);
            D25lag(j)=INMODEL(46);
            D35lag(j)=INMODEL(47);
            D45lag(j)=INMODEL(48);
            A16lag(j)=INMODEL(49);
            A26lag(j)=INMODEL(50);
            A36lag(j)=INMODEL(51);
            A46lag(j)=INMODEL(52);
            D16lag(j)=INMODEL(53);
            D26lag(j)=INMODEL(54);
            D36lag(j)=INMODEL(55);
            D46lag(j)=INMODEL(56);
            sstotApprox1(j)=STATS.SStotal;        %calculate R^2
            ssresidApprox1(j)=STATS.SSresid;
            rsq = 1- ssresidApprox1./sstotApprox1;
            rsq(rsq==Inf) = NaN  ;                   %Set Inf to NaN
            rmse(j)=STATS.rmse;              %Extract rmse
            rmse(rmse==Inf) = NaN;           %Set Inf to NaN
            % repeat regresson only on the sigificant variables
            if sum(INMODEL,2)>0
            xxn = zeros(size(xx,1),8);  
            xip=0; 
                for k=1:8   %8 refers to previous 8 variables including intecept
                    if INMODEL(1,k)==1
                    xip=xip+1;
                    xxn(:,xip)=xx(:,k); %these should now be the correct dimension
        end
    end
              [Bn,SEn,PVALn,INMODELn,STATSn,NEXTSTEPn,HISTORYn]= ...
                                           stepwisefit(xxn,yy,'penter',.05);
                    rmsen(j)=STATSn.rmse;              %Extract rmse
                    rmsen = [rmse,zeros];
                    rmsen(rmse==0) = NaN;           %Set Inf to NaN
            sstotApprox1(j)=STATSn.SStotal;        %calculate R^2
            ssresidApprox1(j)=STATSn.SSresid;
            rsqn = 1- ssresidApprox1./sstotApprox1
            rsqn = [rsqn,zeros];
            rsqn(rsq==0) = NaN 
            yest=xxn*Bn+STATSn.intercept
                R1 = sum((yy-mean(yy)).*(yest-mean(yest)))
                R2 = (sum((yy-mean(yy)).^2)).^.5
                R3 = (sum((yest-mean(yest)).^2)).^.5
                R = (R1/((R2)*(R3)))
                testR(1,j)=R;
            end
          end
        end