% Blasius_Diffop_Example.m
% script solving the nonlinear Blasius equation with discretized
%       differential operators
%
%% inputs
clear all
close all
clc

a = 0;
b = 40;
% mesh points
N = 81;

% Chebyshev grid and differentiation matrix on [-1, 1]
[Dcheb,xcheb] = cheb(N);

% mapping chabyshev grid to the computational domain, non-uniform mapping
yhalf = 6;
c1 = b*yhalf/(b-2*yhalf);
c2 = 1 + 2*c1/(b-a);
dxi = (c2-xcheb).^2 / (c1*(c2+1));
% here are the y coordinates and the differentiation matrix
y = c1*(1+xcheb)./(c2-xcheb) + a;
D1 = diag(dxi) * Dcheb;

%% solving the Blasius equation

[U,dU,ddU] = BlasiusprofDiffop(y,D1);

figure;
plot(U,y,'LineWidth',1.3);
hold on;
plot(dU,y,'LineWidth',1.3,'LineStyle','--');
plot(ddU,y,'LineWidth',1.3,'LineStyle','-.');

xlabel('Velocity');
ylabel('y');
legend({'U','dU','ddU'});

%%
function [U,dU,ddU] = BlasiusprofDiffop(y,D1,varargin)
% function for solving the nonlinear Blasius equation
%       u''' + 0.5*u''*u = 0
% with this, the base flow is U = u', U' = u'', and U'' = u''' = -0.5*u*u''
% inputs:
%       y: y coordinates 
%       D1: differentiation matrix, first derivative
%       additional argument (handles by varargin): initial guess for the solution

% higher derivatives
D2 = D1^2;
D3 = D2*D1;

% indices: one at y=ymax, two at y=0
n = length(y);
rind = [1,n-1,n];
kind = [2:n-2];

% C matrix with boundary conditions
C = [...
    D1(1,:);...              % f'(inf) = 1
    D1(end,:); ...           % f'(0) = 0;
    zeros(1,n-1),1 ...      % f(0) = 0;
    ];

% f: inhomogeneous BC
f = [1;0;0];

% C reordering
C_r = C(:,rind);
C_k = C(:,kind);
% calculating G and H
G = -C_r\C_k;
H = inv(C_r);

% differentiation matrices with boundary conditions - working at points 'k'
D1_k = D1(kind,kind) + D1(kind,rind)*G;
D2_k = D2(kind,kind) + D2(kind,rind)*G;
D3_k = D3(kind,kind) + D3(kind,rind)*G;

% f with boundary conditions - for differentiation - at points 'k'
f_D1_k = D1(kind,rind)*H*f;
f_D2_k = D2(kind,rind)*H*f;
f_D3_k = D3(kind,rind)*H*f;

% function handle for nonlinear equation solver fsolve
fhandle = @(uk) Blassolve(D2_k,D3_k,f_D2_k,f_D3_k,uk);

% initial guess: if not specified, a suitable one is chosen
if nargin < 3
    uk_0 = y(kind);
else
    uk_0 = varargin{1}(kind);
end

% nonlinear equation solver parameters
optim = optimoptions('fsolve','MaxIterations',1000,'FunctionTolerance',1e-12,...
    'StepTolerance',1e-12,'OptimalityTolerance',1e-12,'MaxFunctionEvaluations',20000);

% solution of the Blasius equation
u_k = fsolve(fhandle,uk_0,optim);

% reconstructing the full solution
ufull = zeros(n,1);                 % allocation
u_r = G*u_k + H*f;                  % solution at the 'boundary points'
ufull(kind,1) = u_k;                % reconstruction of the solution using indexing
ufull(rind,1) = u_r;                % reconstruction of the solution using indexing

% Velocity profile obtained by simply differentiation
U = D1*ufull;
dU = D2*ufull;
ddU = D3*ufull;

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rhs = Blassolve(D2_,D3_,f_D2,f_D3,u)
% this function formulates the blasius differential eqaution f''' + f*f'' = 0
% since the boundary conditions are inhomogeneous, the differentiation
% includes and inhomogeneous part
rhs = (D3_*u + f_D3) + 0.5*u.*(D2_*u + f_D2);
end
%% Chebyshev differentiation function
function [D,x] = cheb(N)
if N==0, D=0; x=1; return, end
x = cos(pi*(0:N)/N)';
c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X-X';
D = (c*(1./c)')./(dX+(eye(N+1))); % off-diagonal entries
D = D - diag(sum(D')); % diagonal entries
end