function blasius
% routine solves the Blasius equation
%
disp('Routine solves the Blasius equation')
eta1 = input('Choose the end point eta =');
%
eta0 = 0.0;
[eta,y] = ode45('fblas', [eta0 eta1], [0.0; 0.0; 0.332]);
subplot(2,1,1); plot(eta,y(:,1));
ylabel 'f'
subplot(2,1,2); plot(eta,y(:,2));
xlabel '\eta'
ylabel 'df'

function dy = fblas(t,y)

dy = [y(2); y(3); -y(1)*y(3)/2.0];
 
function blasiu2
%
% routine computes momentum thickness of the Blasius problem
%
global eta0;
global f0;
global df0;
global ddf0;
%
format long;

disp('Routine computes momentum thickness of the Blasius problem')
u1 = input('Choose the end point u =');
%
u0 = 0.0;
theta0 = 0.0;
%
eta0 = 0.0;
eta1 = 1.e-9;
[eta,f] = ode45('fblas', [eta0 eta1], [0.0; 0.0; 0.332]);
%
eta0  = eta1;
f0    = f(size(f,1),1);
df0   = f(size(f,1),2);
ddf0  = f(size(f,1),3);
%
[u,theta] = ode45('fblas2', [u0 u1], [theta0]);
disp('The momentum thickness of Blasius problem at u=5 is');
theta(size(theta,1),1)

function dtheta = fblas2(eta,theta)
global eta0;
global f0;
global df0;
global ddf0;
%
eta1 = eta;
[eta00,f] = ode45('fblas', [eta0 eta1], [f0; df0; ddf0]);
df = f(size(f,1),2);
dtheta = [df*(1.0-df)];
 
