%   Falkner-Skan solution to the boundary layer, computed numerically
%   f'''+f*f''+beta*(1-f'^2)=0  or f'''=-f*f''-beta*(1-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''-beta*(1-f'^2)]
    clear all
    clf reset
    global range
    figure(2)
    clf reset
    beta=[1.5 1 .5  0 -0.198];
    range=linspace(0,8,160);
    soln=[]; betad=[];
    for ibeta=1:length(beta)
    niter=10;
    A=linspace(0,2,niter);    % initial estimate for f''
    flag=0;
    figure(1)
    clf reset
    for iter=1:niter           % compute ODE to verify solution is possible
        clear eta Y
        Vinit=[0 0 A(iter)];
        options=odeset('Events',@events,'RelTol',1.e-6,'AbsTol',1.e-7);
        [eta,Y]=ode45(@falknerf,range,Vinit,options,beta(ibeta));
        plot(eta,Y(:,2),'-')
        hold on
        if flag==0&Y(end,2)>1  % find where f' crosses +1
            iflag=iter;
            flag=1;
        end
    end
    clear eta Y
    figure(2)
    warning off MATLAB:fzero:UndeterminedSyntax
    opt=optimset('Display','iter');
    if iflag==1 
        iflag=2;
    end
    A1=fzero('fpeq1',[A(iflag-1) A(iflag)],[],beta(ibeta));
    Vinit=[0 0 A1];
    [eta,Y]=ode45(@falknerf,range,Vinit,options,beta(ibeta));
    plot(eta,Y(:,2),'-k')
    hold on
    axis([0 6 0 1])
    xlabel('eta')
    ylabel('Dimensionless velocity')
    title(['Falkner-Skan , Beta = ',num2str(beta)])
    soln=[soln eta Y(:,1) Y(:,2)];
    betad=[betad beta(ibeta)];
    end
    save solnd.dat soln -ascii
    save betad.dat beta -ascii
    
    

