X = [ones(length(year), 1), year, year.^2]; % Design matrix with ones, years, and years^2 % Step 2: Create the vector z z = conc; % This is the vector of observed concentrations % Step 3: Create the matrix S S = X' * X; % Normal equation for least squares % Step 4: Find the Cholesky decomposition of S and store it in U U = chol(S, 'lower'); % Perform Cholesky decomposition % Step 5: Solve the triangular system U*y = X'*z y = U \ (X' * z); % First solve for y % Step 6: Solve the triangular system U'*beta = y to get beta beta = U' \ y; % Solve for beta (the coefficients) % Step 7: Compute the best fit values best_fit_values = X * beta; % Calculate polynomial values using the estimated coefficients % Step 8: Plot the fit in red together with the original data points figure; plot(year, conc, 'bo', 'MarkerFaceColor', 'b', 'DisplayName', 'Original Data'); % Original data points hold on; plot(year, best_fit_values, 'r-', 'LineWidth', 2, 'DisplayName', 'Best Fit Quadratic'); % Best fit line hold off; % Finalizing the plot xlabel('Year'); ylabel('CO2 Concentration (ppm)'); title('CO2 Concentration Best Fit Quadratic Polynomial'); legend show; % Show legend axis tight; % Adjust axes tightly around the data grid on; % Optional: Add a grid for better readability