function [u,y]=blasius(L,N); 
% a script to compute the blasius profile, using
% newton solver fsolve, and Chebyshev collocation.
% 	input:
% L: box height
% N: number of gauss-Lobatto points
% 	output:
% u: blasius profile, normalized to deltastar=1
% y: collocation points

%%%% differentiation matrices, scaling from [-1,1] to [0,L]
[yvec,DM] = chebdif(N,3);
scale=-2/L;
D1=DM(:,:,1)*scale;    
D2=DM(:,:,2)*scale^2;    
D3=DM(:,:,3)*scale^3;  	

%%%% prepare boundary conditions
d=[N]; % indices with Dirichlet condition
n=[1 N];% indices with neuman condition, wall is at indice N
r=2:N-2;% remaining degrees of freedom
s=[1 N-1];%removed degree of freedom, due to Neuman
b=[1;0]; % values of neuman conditions
a=[0]; % values of dirichlet conditions

DN=-D1(n,s)\D1(n,r);% dynamic term 
DF=D1(n,s)\(b-D1(n,d)*a); % forcing term 

%%%% implement boundary conditions
DD1=D1(r,r)+D1(r,s)*DN;   D1F=D1(r,d)*a+D1(r,s)*DF;   
DD2=D2(r,r)+D2(r,s)*DN;   D2F=D2(r,d)*a+D2(r,s)*DF;   
DD3=D3(r,r)+D3(r,s)*DN;   D3F=D3(r,d)*a+D3(r,s)*DF;   

%%%% initial guess and solution
f0=-(2*yvec(r)-2)*L/4-1.21.*(1-exp((yvec(r)-1)*L/2));
f=fsolve(@(f)DD3*f+D3F+f.*(DD2*f+D2F),f0,[0,1e-14,1e-14]);

%%%% recover full solution
ff=zeros(N,1);
ff(r)=f; % kept degrees of freedom
ff(s)=DN*f+DF; % removed points
ff(d)=a; % dirichlet points

%%%% integration weigths
nW=0:1:N-1;
jW=0:1:N-1;
bW=ones(1,N);
bW(1)=0.5;
bW(N)=0.5;
cW=2*bW;
bW=bW/(N-1);
S=cos(nW(3:N)'*jW*(pi/(N-1)));
IW=L/2*diag(bW.*[(2+(cW(3:N).*((1+(-1).^nW(3:N))./(1-nW(3:N).^2)))*S)]);
iw=diag(IW);

%%%% normalise to deltastar=1
u=D1*ff;  
deltastar=iw'*(1-u);
y=L*0.5*(yvec+1)/deltastar;
L=L/deltastar;
