function similarity_boundary_layer
    % Initial guess for F''(0)
    beta_guess = 0.5;

    % Define the span of integration
    eta_span = [0, 10];

    % Use fsolve to find the correct initial condition for F''(0)
    options = optimset('Display', 'iter');
    beta = fsolve(@(beta) boundary_conditions(eta_span, beta), beta_guess, options);
    
    % Solve the ODE with the correct initial condition
    [eta, y] = ode45(@ode_system, eta_span, [0, 0, beta]);
    
    % Define parameters for the physical domain
    nu = 1e-5;
    U_inf = 1;
    x = linspace(0.1, 2, length(eta));  % avoid division by zero by starting from 0.1
    y_physical = linspace(0, 10, 100);
    
    % Create a meshgrid for x and y
    [X, Y] = meshgrid(x, y_physical);
    
    % Calculate eta based on x and y
    Eta = Y .* sqrt(U_inf ./ (nu .* X));
    
    % Interpolate F and F' values from the solution of ODE
    F = interp1(eta, y(:, 1), Eta, 'linear', 'extrap');
    F_prime = interp1(eta, y(:, 2), Eta, 'linear', 'extrap');
    
    % Calculate the velocity components u and v
    U = U_inf .* F_prime;
    V = (U_inf / 2) .* (Eta .* F_prime - F);
    
    % Plot streamlines
    figure;
    streamslice(X, Y, U, V);
    xlabel('x');
    ylabel('y');
    title('Streamlines');
    
    
    % Plot the results of F, F', F''
    figure;
    plot(eta, y(:,1), '-r', 'DisplayName', 'F(\eta)');
    hold on;
    plot(eta, y(:,2), '-g', 'DisplayName', 'F''(\eta)');
    plot(eta, y(:,3), '-b', 'DisplayName', 'F''''(\eta)');
    legend;
    title('Solution of 2F'''''' + FF'''' = 0 using the Shooting Method');
    xlabel('\eta');
    ylabel('Solution');
    grid on;
        % Parameters
    Pr = 1; % Prandtl number
    f =0.3; % Example function f, you should replace this with your specific function
    
    % Shooting method parameters
    eta_span = [0, 10]; % Integration span for eta
    theta_guess = 0.5; % Initial guess for theta'(0)
    
    % Use fsolve to find the correct initial condition for theta'(0)
    theta_prime_0 = fsolve(@(theta_prime_0) boundary_conditionsH(eta_span, theta_prime_0, Pr, f), theta_guess);
    
    % Solve the ODE with the correct initial condition
    [eta, sol] = ode45(@(eta, y) ode_systemH(eta, y, Pr, f), eta_span, [0; theta_prime_0]);
    theta = sol(:,1);
    theta_prime = sol(:,2);
    theta_double_prime = -Pr/2 * f * theta_prime;
    
    % Plot results
    figure;
    plot(eta, theta_double_prime, '-r', 'DisplayName', '\theta''''(\eta)')
    hold on
    plot(eta, theta_prime, '-g', 'DisplayName', '\theta''(\eta)');
    hold on
    plot(eta, theta, '-b', 'DisplayName', '\theta(\eta)');
    hold on
    xlabel('\eta');
    ylabel('\theta(\eta)');
    legend;
    grid on;
    figure;
    plot(x,3*y(:,2),'linewidth',2)
    xlabel('x(m)')
    ylabel('u (m/s)')
    title('u  vs. x')
    grid on
    theta = exp(-eta); % Temperature profile in similarity variable
    T = interp1(eta, theta, Eta(:), 'linear', 'extrap'); % Interpolated temperature field
    
    % Reshape T back into the grid shape
    T = reshape(T, size(Eta));

    % Calculate heat function H (Paulhasen method)
    % H = integral(U*dT - V*dX)
    H = zeros(size(X));
    for i = 2:size(X, 1)
        for j = 2:size(X, 2)
            H(i, j) = H(i-1, j) + U(i, j) * (T(i, j) - T(i-1, j)) * (y_physical(2) - y_physical(1)) - ...
                      V(i, j) * (X(i, j) - X(i, j-1)) * (x(2) - x(1));
        end
    end

    % Plot heatlines using contour
    figure;
    contour(X, Y, H, 20); % 20 contour levels
    xlabel('x');
    ylabel('y');
    title('Heatlines');
    
    
end

function dydeta = ode_system(~, y)
    % Define the system of first-order ODEs
    % y(1) = F, y(2) = F', y(3) = F''
    dydeta = zeros(3,1);
    dydeta(1) = y(2);
    dydeta(2) = y(3);
    dydeta(3) = -0.5 * y(1) * y(3);
end

function res = boundary_conditions(eta_span, beta)
    % Solve the ODE with the initial guess for F''(0)
    [~, y] = ode45(@ode_system, eta_span, [0, 0, beta]);
    
    % The boundary condition we want to satisfy is F'(\eta_end) = 1
    eta_end = eta_span(end);
    F_prime_at_eta_end = y(end, 2);
    res = F_prime_at_eta_end - 1;
end
function dydeta = ode_systemH(eta, y, Pr, f)
    % y(1) = theta, y(2) = theta'
    dydeta = zeros(2,1);
    dydeta(1) = y(2);
    dydeta(2) = -Pr/2 * f * y(2);
end

function res = boundary_conditionsH(eta_span, theta_prime_0, Pr, f)
    % Solve the ODE with the initial guess for theta'(0)
    [~, sol] = ode45(@(eta, y) ode_systemH(eta, y, Pr, f), eta_span, [0; theta_prime_0]);
    
    % The boundary condition we want to satisfy is theta'(\infty) = 1
    theta_prime_at_infinity = sol(end, 2);
    res = theta_prime_at_infinity - 1;
end