%%%% parameters
N=100; % number of points in y
Re=1000; % Reynolds number based on deltastar
alpha=0; % streamwise wavenumber
beta=0.7; %spanwise wavenumber
L=15; % box height 

%%%% compute blasius profile and interpolate
[ubla,ybla]=blasius(20,150); % compute Blasius profile
if ybla(1)<L; disp(['Error... increase L']);break;end;
yto =-(1-(L/2*(1-cos((N+1:-1:0)*pi/(N+1))))*2/ybla(1));
u=chebint(ubla,yto); % interpolate to new grid

%%%% compute OSS
[A,B,C,Q,yvec]=oss(N,alpha,beta,Re,L,u);

%%%% Solve differential Lyapunov equation
f=[zeros(N,1);exp(-((yvec-2)/0.7).^2);zeros(N,1)];% forcing function
M=f*f';% covariance of random forcing
P0=zeros(2*N,2*N);% covariance of initial condition
tmax=5000;nt=20;% time vector
[PP,rms,tvec]=difflyap(A,B,C,M,P0,Q,nt,tmax);% time evolution of P
mesh(tvec,yvec,rms(1:N));%plot rms

