function [A,B,C,Q,y]=oss(N,alpha,beta,Re,L,u)
% build the Orr-Sommerfeld/Squire system using Chebyshev collocation
%        inputs:
% N: number of inner points
% alpha, beta: streamwise and spanwise wavenumbers
% Re: Reynolds number, based on deltastar
% L,u: boxheight and base flow profile
%        outputs:
% A,B,C: state space operators
% Q, y: inner product matrix and collocation points

%%%% differentiation matrices and base flow
scale=2/L;
[y,DM] = chebdif(N+2,2); %Chebyshev collocation
DM(:,:,1)=DM(:,:,1)*scale;    
DM(:,:,2)=DM(:,:,2)*scale^2;    
y=(y(2:end-1)+1)/scale; 
up=DM(:,:,1)*u; % differentiate base flow
upp=DM(:,:,2)*u;

%%%% implement homogeneous boundary conditions 
sel=2:N+1;
D1=DM(sel,sel,1);
D2=DM(sel,sel,2);
[y,D4]=cheb4c(N+2); % fourth order differentiation matrix
D4=D4*(2/L)^4;

%%%% laplacian
I=eye(N);Z=zeros(N,N);
k2=alpha^2+beta^2;
delta=(D2-k2*I); % Laplacian
delta2=(D4-2*k2*D2+k2*k2*I); % square of the Laplacian

%%%% build oss matrix
LOS=-i*alpha*diag(u(sel))*delta+i*alpha*diag(upp(sel))+delta2/Re;     
LSQ=-i*alpha*diag(u(sel))+delta/Re;
LC= -i*beta*diag(up(sel));
A=[delta\LOS,Z;LC,LSQ]; % the dynamic matrix for OSS

%%%% input,  output, and inner product  operators
B=[delta\(-i*alpha*D1),delta\(-k2*I),delta\(-i*beta*D1);i*beta*I,Z,-i*alpha*I];
C=[i*alpha*D1/k2,-i*beta*I/k2 ; I,zeros(N,N) ; i*beta*D1/k2,i*alpha*I/k2];
Q=enermat(N,L,DM,k2);
