Press "Enter" to skip to content

Solving the Blasius boundary layer equation using Matlab

Boundary layers are present any time a viscous fluid (such as air) is in contact with a surface and there is relative motion between them. Perhaps the simplest example of a boundary layer is flow over a flat-plate at 0 angle of attack. For flat plate flow, you can start with the full Navier-Stokes equations and simplify them to a 3rd-order nonlinear ODE. This is often referred to as the Blasius equation (and its solution as the Blasius boundary layer) after the guy who first discovered it:
$$f^{\prime\prime\prime} + \frac{1}{2}ff^{\prime\prime} = 0.$$
The boundary conditions are no penetration:
$$f(0) = 0,$$
no slip:
$$f^\prime(0) = 0,$$
and the far field must match the freestream flow:
$$\lim_{\eta\to\infty} \:f^\prime(\eta) = 1.$$

The purpose of this post is not to detail the derivation but to show how you can obtain solutions to this ODE using a shooting method and Matlab.

It is the last boundary condition that makes solving this equation a little difficult. We would like to start at \(\eta=0\) and integrate by marching towards the freestream. However to do this we must chose a value of \(f^{\prime\prime}(0)\). Unless we choose this value correctly the far-field boundary condition will not be satisfied.

One solution to this is to make a guess for \(f^{\prime\prime}(0)\), solve the ODE, check the boundary condition, and repeat using the bisection method until the boundary condition is satisfied to within a given tolerance. This is the method implemented in the script below:

function [eta, f] = blasius(eta_max)
  % Initial guesses for f''(0), these values
  % need to bracket the true value
  a = 0.3;
  b = 1.0;
  % The boundary layer only extends to around
  % eta = 5, we don't want to integrate more
  % than necessary (also this sometimes
  % becomes unstable if this value is set too high).
  % Also note that this is the height where the
  % freestream boundary condition is enforced.
  eta_max_ode = 20;
  % Maximum iterations
  max_iter = 100;
  % Iteration counter
  count = 0;
  % Tolerance for boundary condition
  epsilon = 1e-10;
  % Value greater than epsilon to start it off
  error = 1e15; 
  while(abs(error) > epsilon && count < max_iter)
    % Solve the ODE using ode45
    sigma = (a+b)/2;
    f0 = [0, 0, sigma];
    [eta_ode, f_ode] = ode45(@BlasiusFunc, [0, eta_max_ode], f0);
    % Perform bisection method
    error = f_ode(end,2) - 1;
    if(error > 0)
      b = sigma;
      a = sigma;
    % Increment the counter
    count = count + 1;
  f = f_ode;
  eta = eta_ode;
  % Adjust the ends of f and eta to respect eta_max
  if(eta_max > eta_max_ode)
    f = [f; [(eta_max - eta_max_ode + f(end,1)), 1, 0]];
    eta = [eta; eta_max];
  elseif(eta_max < eta_max_ode)
    f = f(eta < eta_max,:);
    eta = eta(eta < eta_max);
    eta = [eta; eta_max];
    f_interp = interp1(eta_ode, f_ode, eta_max);
    f = [f; f_interp];

function [df] = BlasiusFunc(eta, f)
  df = zeros(3,1);
  df(1) = f(2);
  df(2) = f(3);
  df(3) = -f(1)*df(2)/2;

And here are some plots of the solutions:
Blasius solution for f.
Blasius solution for f‘.
Blasius solution for f‘‘.


  1. Ferdousi Ema
    Ferdousi Ema 2017-08-11

    Hello brother,
    I used to run the code that u gave above
    but it showed a error
    “” ??? Input argument “eta_max” is undefined.

    Error in ==> blasius2 at 48
    if(eta_max > eta_max_ode)
    now what can i do ?
    please tell me about the problem

    • jason
      jason 2018-09-25

      The code above is just a function. To use it save it to a file called blasius.m and call it from either the command window or another script as:
      [eta, f] = blasius(10);

      You can then plot it with:
      plot(f(:,2), eta);
      xlabel('$$\frac{u}{U_\infty}$$', 'interpreter', 'latex');
      ylabel('$$\eta$$', 'interpreter', 'latex');

    SHIBU CLEMENT 2017-10-24

    Thanks Jason for posting the Matlab code for solving Blasius equation using Runge-Kutta and Bisection methods. Instead of bisection method if we use Newton-raphson method, what changes we have to make it in the code?


    • jason
      jason 2018-09-25

      Interesting idea! The biggest challenge will be to find the derivative of R = f'(\eta_{max}) - 1 with respect to f''(0). Once you find that you can apply simple 1D Newton method to find the f''(0) that makes R zero.

  3. Golam Mortuja
    Golam Mortuja 2018-08-14

    the program is not running in MatLab 2017a. Is the erro=1e15 correct(in line 25)?

    • jason
      jason 2018-09-25

      Hi Golam, I just ran the code in MATLAB 2017b without errors. What problem are you seeing?

Leave a Reply

Your email address will not be published. Required fields are marked *