inflectionPointTestPoly

%Testing code that calculates inflection points for a polynomial function

% Define x and y for either a known function or numerical data
% Example for a known function y = x^5 + 5*x^3
x = linspace(-10, 10, 5000); % Higher resolution for better detection
y = x.^7 - 7*x.^5 + 14*x.^3 - 7*x; % Define the function

% 1. Interpolate and smooth the data for higher resolution
x_interp = linspace(min(x), max(x), 10000); % High-resolution x values
y_interp = interp1(x, y, x_interp, 'pchip'); % Smooth interpolation
y_smooth = smoothdata(y_interp, 'gaussian', 50); % Stronger Gaussian smoothing

% 2. Calculate first and second derivatives
dy_dx = gradient(y_smooth) ./ gradient(x_interp); % First derivative
d2y_dx2 = gradient(dy_dx) ./ gradient(x_interp); % Second derivative

% 3. Detect zero-crossings in the second derivative
sign_changes = find(diff(sign(d2y_dx2)) ~= 0); % Detect strict sign changes

% 4. Initialize arrays for inflection points, ensuring unique detections
inflection_points_x = [];
inflection_points_y = [];

for i = 1:length(sign_changes)
    inflection_x = x_interp(sign_changes(i));
    inflection_y = y_smooth(sign_changes(i));
    
    % Ignore points near the edges, e.g., within 5% of each end of the x range
    if inflection_x > x_interp(round(0.05 * length(x_interp))) && ...
       inflection_x < x_interp(round(0.95 * length(x_interp)))
    
        % Ensure unique inflection points based on proximity
        if isempty(inflection_points_x) || all(abs(inflection_x - inflection_points_x) > 1e-3)
            inflection_points_x = [inflection_points_x, inflection_x];
            inflection_points_y = [inflection_points_y, inflection_y];
            
            % Display each unique inflection point
            disp(['Inflection point: (x, y) = (' num2str(inflection_x) ', ' num2str(inflection_y) ')']);
        end
    end
end

% 5. If no inflection points are found, display a message
if isempty(inflection_points_x)
    disp('No inflection points found.');
end