%   Blasius solution to the boundary layer, computed numerically
%   2f'''+f*f''=0  or f'''=-f*f''/2
%   f is function of similarity variable, eta
%   f(0)=0, f'(0)=0, f'(infinity)=1
%   f''(0)= to be estimated to satisfy condition at infinity
%   Y=[f f' f''],  dY/deta=[f' f'' f''']=[Y(2) Y(3) -f*f''/2]
    clear all
    load blasiusd.dat;
    A=1;    % initial estimate for f''
    Vinit=[0 0 A];
    range=[0 8];
    for iter=1:2
        if iter>1
            A=A*Y(end,2)^(-3/2);
            Vinit=[0 0 A];
        end
        [eta,Y]=ode45('blasiusf',range,Vinit);
        figure(iter)
        clf reset
        plot(eta,Y(:,2),'-k',blasiusd(:,1),blasiusd(:,2),'g*')
        xlabel('eta')
        ylabel('Dimensionless velocity')
        title('Numerical Solution of Blasius Equation')
        Y(end,2)
    end
