%clear all; close all;

global beta N

%% grid
%N = 180; % degree of highest Chebyshev polynomial
%[D,y]=cheb(N); %from Trefethen, Spectral Methods in Matlab
%D2 = D*D; 
%S = diag([0; 1./(1-y(2:N).^2); 0]);
%D4 = (diag(1-y.^2)*D^4 - 8*diag(y)*D^3 - 12*D^2)*S;

[D,y]=cheb(N); %from Trefethen, Spectral Methods in Matlab
D2 = D*D; 
S = diag([0; 1./(1-y(2:N).^2); 0]);
D4 = (diag(1-y.^2)*D^4 - 8*diag(y)*D^3 - 12*D^2)*S;
%D2=D2(2:N,2:N);
%D4=D4(2:N,2:N);
%D=D(2:N,2:N);

H = 10; 
y = (y+1)*H; 
D = D/H; 
D2=D2/(H^2); 
D4=D4/(H^4);
%y=y+1;
%y = y(2:N);

% mean flow
% Blasius boundary layer
delta = 1.;
fact = 1./delta * 5.;

beta = 0; % Falkner-Skan; 
% put beta=0 for Blasius
% put beta>0 for favourable pressure gradient
% put beta<0 for adverse pressure gradient

options = optimset('TolX',1e-10);
ypp = fminsearch(@findfpp,0.4,options);

f0 = [0 0 ypp];
[xout,yout]=ode45(@blasius,y(N+1:-1:1),f0);
%[xout,yout]=ode45(@blasius,y(N-1:-1:1),f0);

U = yout(N+1:-1:1,2);
%U = yout(N-1:-1:1,2);


%% find boundary layer thickness and normalise y
delta = interp1(U,y,0.99);
factor = 1/delta;
y = y*factor;

D = (1/factor)*D;
D2 = (1/factor^2)*D2;
D4 = (1/factor^4)*D4;

DU=D*U;
D2U=D2*U;

figure;
plot(U,y,'b');
hold on
plot(DU,y,'r');
plot(D2U,y,'k');

xlabel('U');ylabel('y/\delta'); grid on;
%y=y(2:N);
%U=U(2:N);