Chebyshev polynomials in Matlab

Hi. I am now trying to solve stability problem in rotating disk flow. I have to use Chebyshev spectral methods in Matlab in order to get the solutions and graphs. But, i have one problem about this method. The details are as below:

Basically, what I need to do and what I have done:

1) Solve the governing equations and transform into non-linear ODE (done)

2) Solve the non-linear ODE from eta=0 to eta=20

2*U+(1-n/n+1)*eta*U'+W'=0,
U^2-(V+1)^2+(W+(1-n/n+1)*eta*U)*U'-(nu*U')'=0,
2*U*(V+1)+(W+(1-n/n+1)*eta*U)*V'-(nu*V')'=0,
boundary conditions: U(0)=V(0)=W(0)=0, U(infinity)=0, V(infinity)=-1 where infinity=20
n=power law index=from 0.6 to 1
nu=(U'^2+V'^2)^(n-1)/2

I solve above non-linear equations by maple and I got th numerical values for U(0) until U(20), same goes to V and W.

3) I have to perturb the equations thentransform those equations into normal modes form then do some linearisation (done)

4) HERE THE PROBLEM ARISE. In order to do chebyshev spectral methods in Matlab, i have to transform the linear disturbance equations into Chebyshev polynomials form (which I have successfully done). The problem is, in the transformed equations, there is U, V and W which is the mean velocity where I dont know how to assign these functions. Below is the code:

function [U,V,W,P,C,y,C_ls]=orr_somm_rotating_channel_simple(alpha,beta,Re,omega)
  N=6; % number of collocation points
  n=0:N;
  y=cos(pi/N*n); % gauss-Lobatto points

    % The recurrence relations of the trigo func
    T(1,1:(N+1))=1;
    T(2,1:(N+1))=y;
    for k=3:N+1
        for ll=1:N+1
            T(k,ll)=2*y(ll)*T(k-1,ll)-T(k-2,ll);
        end
    end

    % The derivatives of Chebyshev polynomials
    T1(1,1:(N+1))=0;
    T1(2,:)=T(1,:);
    T1(3,:)=4*T(2,:);
    for k=3:N
        for ll=1:N+1             
            T1(k+1,ll)=2*1*T(k,ll)+2*y(ll)*T1(k,ll)-T1(k-1,ll);
        end
    end

    T2(1,1:(N+1))=0;
    T2(2,:)=T1(1,:);
    T2(3,:)=4*T1(2,:);
    for k=3:N
        for ll=1:N+1
            T2(k+1,ll)=2*2*T1(k,ll)+2*y(ll)*T2(k,ll)-T2(k-1,ll);
        end
    end

    T=T';
    T1=T1';
    T2=T2';

    % flow variables
     dydr=1;
         i=sqrt(-1);
     er=-200*i;

       for n=0:N
          for k=1:N-1 % I skip only the r=1 point. The BC will be posed at the end of the program
          A(4*k+1,4*n+1)=(U(k+1)*i*alpha+1/Re*(alpha^2+beta^2))*T(k+1,n+1)-1/Re*T2(k+1,n+1)*dydr^2;
          A(4*k+1,4*n+2)=(U1(k+1)-omega)*T(k+1,n+1);
          A(4*k+1,4*n+3)=0;
          A(4*k+1,4*n+4)=i*alpha*T(k+1,n+1);

          A(4*k+2,4*n+1)=omega*T(k+1,n+1);
          A(4*k+2,4*n+2)=(V(k+1)*i*alpha+1/Re*(alpha^2+beta^2))*T(k+1,n+1)-1/Re*T2(k+1,n+1)*dydr^2;
          A(4*k+2,4*n+3)=0;
          A(4*k+2,4*n+4)=T1(k+1,n+1)*dydr;

          A(4*k+3,4*n+1)=0;
          A(4*k+3,4*n+2)=0;
          A(4*k+3,4*n+3)=(W(k+1)*i*alpha+1/Re*(alpha^2+beta^2))*T(k+1,n+1)-1/Re*T2(k+1,n+1)*dydr^2;
          A(4*k+3,4*n+4)=i*beta*T(k+1,n+1);

          A(4*k+4,4*n+1)=i*alpha*T(k+1,n+1);
          A(4*k+4,4*n+2)=T1(k+1,n+1)*dydr;
          A(4*k+4,4*n+3)=i*beta*T(k+1,n+1);
          A(4*k+4,4*n+4)=0;

          B(4*k+1,4*n+1)=i*T(k+1,n+1);
          B(4*k+1,4*n+2)=0;
          B(4*k+1,4*n+3)=0;
          B(4*k+1,4*n+4)=0;

          B(4*k+2,4*n+1)=0;
          B(4*k+2,4*n+2)=i*T(k+1,n+1);
          B(4*k+2,4*n+3)=0;
          B(4*k+2,4*n+4)=0;

          B(4*k+3,4*n+1)=0;
          B(4*k+3,4*n+2)=0;
          B(4*k+3,4*n+3)=i*T(k+1,n+1);
          B(4*k+3,4*n+4)=0;

          B(4*k+4,4*n+1)=er*i*alpha*T(k+1,n+1);
          B(4*k+4,4*n+2)=er*T1(k+1,n+1)*dydr;
          B(4*k+4,4*n+3)=er*i*beta*T(k+1,n+1);
          B(4*k+4,4*n+4)=0; 

          end

A(1,4*n+1)=er*T(1,n+1); A(1,4*n+2)=0; A(1,4*n+3)=0; A(1,4*n+4)=0;

A(2,4*n+1)=0; A(2,4*n+2)=er*T(1,n+1); A(2,4*n+3)=0; A(2,4*n+4)=0;

          A(3,4*n+1)=0;                 A(3,4*n+2)=0;                    A(3,4*n+3)=er*T(1,n+1);        A(3,4*n+4)=0;
       
          A(4,4*n+1)=0;                 A(4,4*n+2)=er*T1(1,n+1);         A(4,4*n+3)=0;                  A(4,4*n+4)=0;

          B(1,4*n+1)=T(1,n+1);       B(1,4*n+2)=0;                B(1,4*n+3)=0;               B(1,4*n+4)=0;  

         B(2,4*n+1)=0;              B(2,4*n+2)=T(1,n+1);         B(2,4*n+3)=0;               B(2,4*n+4)=0;

          B(3,4*n+1)=0;              B(3,4*n+2)=0;                B(3,4*n+3)=T(1,n+1);        B(3,4*n+4)=0;

          B(4,4*n+1)=0;              B(4,4*n+2)=T1(1,n+1);        B(4,4*n+3)=0;               B(4,4*n+4)=0;

          A(4*N+1,4*n+1)=er*T(N+1,n+1);       A(4*N+1,4*n+2)=0;                      A(4*N+1,4*n+3)=0;                      A(4*N+1,4*n+4)=0;

A(4*N+2,4*n+1)=0; A(4*N+2,4*n+2)=er*T(N+1,n+1); A(4*N+2,4*n+3)=0; A(4*N+2,4*n+4)=0;

          A(4*N+3,4*n+1)=0;                   A(4*N+3,4*n+2)=0;                      A(4*N+3,4*n+3)=er*T(N+1,n+1);          A(4*N+3,4*n+4)=0;
       
          A(4*N+4,4*n+1)=0;                   A(4*N+4,4*n+2)=er*T1(N+1,n+1);         A(4*N+4,4*n+3)=0;                      A(4*N+4,4*n+4)=0;

          B(4*N+1,4*n+1)=T(N+1,n+1);          B(4*N+1,4*n+2)=0;                   B(4*N+1,4*n+3)=0;                   B(4*N+1,4*n+4)=0;  

         B(4*N+2,4*n+1)=0;                   B(4*N+2,4*n+2)=T(N+1,n+1);          B(4*N+2,4*n+3)=0;                   B(4*N+2,4*n+4)=0;

          B(4*N+3,4*n+1)=0;                   B(4*N+3,4*n+2)=0;                   B(4*N+3,4*n+3)=T(N+1,n+1);          B(4*N+3,4*n+4)=0;

          B(4*N+4,4*n+1)=0;                   B(4*N+4,4*n+2)=T1(N+1,n+1);         B(4*N+4,4*n+3)=0;                   B(4*N+4,4*n+4)=0;
    end

    C1=zeros(4*(N+1),1);
    V1=zeros(4*(N+1),4*(N+1));
    U2=zeros((N+1),4*(N+1));
    V2=zeros((N+1),4*(N+1));
    W2=zeros((N+1),4*(N+1));
    P2=zeros((N+1),4*(N+1));
    [V,C]=eig(A,B);
    for kk=1:4*(N+1)
       C1(kk)=C(kk,kk);
    end
    clear C
     C=C1;
    [pippo,paperino]=sort(-imag(C1));
    C=C1(paperino);
    for kk=1:length(paperino)
       V1(:,kk)=V(:,paperino(kk));
    end

% % % % for kk=1:length(paperino) for jj=1:(N+1) U2(:,kk)=U2(:,kk)+V1(1+(jj-1)*4,kk)*T(:,jj); V2(:,kk)=V2(:,kk)+V1(2+(jj-1)*4,kk)*T(:,jj); W2(:,kk)=W2(:,kk)+V1(3+(jj-1)*4,kk)*T(:,jj); P2(:,kk)=P2(:,kk)+V1(4+(jj-1)*4,kk)*T(:,jj); end end U=zeros(N+1,0); V=U; W=U; P=U; %cleaning the eigenmatrices k=1; C1=0; qqq=sum(abs(U2))+sum(abs(V2))+sum(abs(W2)); for j=1:length(paperino) if ((abs(C(j))<2)&(qqq(j)>0))%cut off filter omega<2

          U(:,k)=U2(:,j); 
          V(:,k)=V2(:,j);
          W(:,k)=W2(:,j);
          P(:,k)=P2(:,j);
          C1(k,1)=C(j,1);
           if  (imag(C(j))<(imag(1/er)*(1+1e-3)) && imag(C(j))>(imag(1/er)*(1-1e-3)))
              k=k;
          else
            k=k+1;
          end
       end

    end
    C=C1;
    C_ls=C(1);