% Blasius Equation Solver using Shooting Method

% Define the Blasius ODE as a system of first-order ODEs
function dFdx = blasius_ode(x, F)
    dFdx = [F(2);
            F(3);
            -0.5*F(1)*F(3)];
end

% Define initial guess for shooting method
function eta = blasius_shooting(beta)
    % Initial conditions
    F0 = [0; 0; beta]; % F(0) = 0, F'(0) = 0, F''(0) = beta
    xspan = [0 10]; % From eta=0 to eta=10 (an arbitrary large value)
    
    % Solve the ODE using MATLAB's ODE45 solver
    [x, F] = ode45(@blasius_ode, xspan, F0);
    
    % The correct beta will give F'(eta_inf) = 1
    eta = F(end, 2) - 1;
end

% Main function to find the solution
function blasius_solver
    % Use fzero to find the correct value of beta (F''(0))
    beta_correct = fzero(@blasius_shooting, [0.1, 1.0]); % Initial guess range
    
    % Solve the Blasius ODE using the correct beta
    xspan = [0 10];
    F0 = [0; 0; beta_correct];
    [x, F] = ode45(@blasius_ode, xspan, F0);
    
    % Plot the results
    figure;
    plot(x, F(:, 2), 'b', 'LineWidth', 2);
    xlabel('\eta');
    ylabel('f''(\eta)');
    title('Blasius Solution for Flow over a Flat Plate');
    grid on;
    
    % Display final value of F''(0)
    disp(['The value of F''''(0) is: ', num2str(beta_correct)]);
end

% Run the main function
blasius_solver();