function mnhf_blasius(beta0)
%MNHF_BLASIUS solves the Blasius boundary layer equation using a shooting 
% method. 
%
% beta0 - Guess for the value of \beta(\eta=0) that satisfies
% \Omega(\eta=\infty) \to 1
%
% Usage  mnhf_blasius(0.332)


%% Specify maximum value for eta, the independent variable.
eta_max = 1e1;

%% Specify initial conditions then solve ODE using ode45.
f(1) = 0.0;
Omega(1) = 0.0;
beta(1) = beta0;
[eta,y] = ode45(@blasius_func,[0.0 eta_max],[f(1) Omega(1) beta(1)]);

%% Plot/print results.
figure(1); hold on; box on; set(gca,'fontsize',18)
plot(eta,y(:,1),'k-','linewidth',2)
plot(eta,y(:,2),'k--','linewidth',2)
ylim([0.0 2.0])
xlabel('{\it \eta}'); ylabel('{\it f} (solid), {\it \Omega} (dashed)')
fprintf('Terminal value of Omega: %12.8f\n',y(end,2))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dy=blasius_func(~,y)
%% Specify the three first order equations that are equivalent to the third
%% order Blasius boundary layer equation.

dy = zeros(3,1); % Initialize.
dy(1) = y(2);
dy(2) = y(3);
dy(3) = -0.5*y(1)*y(3);