%{
  Solution of 3rd-oder BVP describing 
  the laminar velocity boundary layer
  over a flat plate.

  Method: shooting, combined with RK4 for
  integration and secant method for zero-search.

  M. Piller, 21st march 2023
%}
function Blasius

% Function handle to system of ODEs (Blasius' problem)
f    = @BlasiusFunc;
% Step size for integration
h0   = 0.05;
% Initial value of independent variable
% (x0 = 0 means "at the wall")
x0   = 0;
% "Numerical" boundary layer extends up to x = 10
xend = 10;
% Tollerance for zero-fnding procedure
toll = 1e-6;

% First two guesses for the second derivative of
% the solution at the wall (needed to start the 
% secant method). We look for the value of the
% second derivative at the wall that yields 
% a value of 1 for the first derivative at 
% x = xend (10 in this case)
A = 0;
B = 1;

[~,y]=RK4stage(x0,xend,[0;0;A],f,h0);
eA = y(2)-1;

[~,y]=RK4stage(x0,xend,[0;0;B],f,h0);
eB = y(2)-1;

% Switch first guesses such that the best one is "B"
if (abs(eA)<abs(eB))
    tmp =  B;  B =  A;  A = tmp;
    tmp = eB; eB = eA; eA = tmp;
end

% Store values
Cstore = [B,eB+1];

% Start loop for secant method
while 1

    % New iterate
    C = A-eA*(B-A)/(eB-eA);
    [~,y]=RK4stage(x0,xend,[0;0;C],f,h0);
    eC = y(2)-1;
    A = B; eA = eB; B = C; eB = eC;

    Cstore = [Cstore;[C,eC+1]];

    % Exit condition
    if (abs(eC)<toll)
        break;
    end
end

% Plot dependence of first derivative at xend
% on second derivative at x0 (progress of the 
% secant method)
figure
plot(Cstore(:,1),Cstore(:,2),'ro-');
grid on

% Plot solutioon of Blasius problem
x = x0:(xend-x0)/50:xend;
y = zeros(3,length(x));
y(:,1) = [0;0;C];
for j=2:length(x)
  [~,y(:,j)]=RK4stage(x(j-1),x(j),y(:,j-1),f,h0);
end

figure
plot(x,y(2,:),'b-');
grid on
box on
xl = xlabel('y/\delta');
yl = ylabel('u/u_\infty');

end

% This is the integration from x0 to xend with RK4
function [x,y]=RK4stage(x,xend,y,f,h)

while 1
    if (x+h)>xend
        h=xend-x;
    end
    [x,y]=RK4_1Step(x,y,f,h,{});
    if x>=xend
        break;
    end
end


end