% Blasius_Shooting_Example.m
clear all
close all
clc

%%

% computational domain
domain = [0, 40];
% 
%% initial guess for the nonlinear problem
input0 = 10;

%% function handle for the nonlinear equation solver
ftosolve = @(input) BlasiusSolve(@(y,u)BlasiusODE(y,u),domain,input);

%% solution of the nonlinear problem - the unknown first derivative
input_sol = fsolve(ftosolve,input0);

%% calculating the solution with the correct initial condition

% ode option object for ODE tolerances
odeopt = odeset('Reltol',1e-12,'Abstol',1e-12);

% initial condition for the ODE
Icond = [0; 0; input_sol];

[y,u_sol] = ode45(@(y,u)BlasiusODE(y,u),domain,Icond,odeopt);

% U = u', U' = u'', and U'' = u''' = -0.5*u*u''
U = u_sol(:,2);
dU = u_sol(:,3);
ddU = -0.5*u_sol(:,1).*u_sol(:,3);

%% plotting
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');
axis([ -0.2, 1.2, 0, 40.]);
legend({'U','dU','ddU'});


%%
function out =  BlasiusSolve(fhandle,domain,ddu0)
% function for solving the nonlinear Blasius equation by the shooting method
%       u''' + 0.5*u''*u = 0
% with this, the base flow is U = u', U' = u'', and U'' = u''' = -0.5*u*u''

% ode option object for ODE tolerances
odeopt = odeset('Reltol',1e-12,'Abstol',1e-12);

% initial condition for the ODE
Icond = [0; 0; ddu0];

[~,sol] = ode45(fhandle,domain,Icond,odeopt);

% derivative of u must be 1 at the far field
out = sol(end,2)-1;


end

function dfdt = BlasiusODE(y,u)
% ode for the Blasius equation
% u1 = u
% u2 = u'
% u3 = u''
% the ODE is
% u1' = u2
% u2' = u3
% u3' = -0.5*u1*u3
dfdt = [ u(2,1); ...
         u(3,1); ...
         -0.5*u(3,1)*u(1,1)];
end