Problem in finding the correct Chebyshev interpolation line

Hi, I have problem in finding the correct Chebyshev interpolation line. I have inputted the calculation procedure in accordance with the Chebyshev theorem. However the regression line is apparently does not agree with the trend of the provided data point. Please let me know I missed something or how to improve my code. Your help is highly appreciated.

%This program is for governing a polinomial equation using Chebyshev
%interpolation formula
clc;  clear all; close all;
t=tic;
%% Data input
x_data = [1830; 1920; 1940; 1975; 2000; 2020]; %x values
y_data = [0.089; 3.52; 4.86; 17.05; 25.5; 35.01]; %y values

%% Chebyshev interpolating formula
% Normalize x-values
x_min = min(x_data);
x_max = max(x_data);
x_data_norm = 2 * (x_data - x_min) / (x_max - x_min) - 1;

% Number of data points
n = length(x_data);

% Compute Chebyshev nodes
x_nodes = cos((2*(1:n) - 1) * pi / (2*n));

% Compute Chebyshev coefficients
a = zeros(n, 1);
for i = 1:n
    T = cos((i-1) * acos(x_nodes));
    a(i) = (2/n) * sum(y_data .* T');
end

% Define the Chebyshev polynomial function
T = @(k, x) cos((k-1) * acos(x));

% Evaluate the interpolating polynomial
x_query = linspace(-1, 1, 100); % Query points
y_interp = zeros(size(x_query));
for i = 1:n
    y_interp = y_interp + a(i) * T(i, x_query);
end

%% Plotting the data and show polynomial equation
% Plot the given data points and the interpolating polynomial
% Construct the interpolation equation
interpolation_eqn = '';
for i = 1:n
    if i == 1
        interpolation_eqn = [interpolation_eqn num2str(a(i))];
    else
        interpolation_eqn = [interpolation_eqn ' + ' num2str(a(i)) ' * T' num2str(i-1) '(x)'];
    end
end

disp(['Interpolation Equation:']);

Interpolation Equation:

disp(['P(x) = ' interpolation_eqn]);

P(x) = 28.6763 + -17.4761 * T1(x) + 3.8073 * T2(x) + -0.17701 * T3(x) + -0.17183 * T4(x) + -1.7569 * T5(x)

% Construct the interpolation equation in terms of powers of x
expanded_eqn = '';
for i = 1:n
    if i == 1
        expanded_eqn = [expanded_eqn num2str(a(i))];
    else
        expanded_eqn = [expanded_eqn ' + ' num2str(a(i)) ' * (2 * cos(' num2str(i-1) '* acos(x)) - 1)^' num2str(i-1)];
    end
end

disp(['Expanded Interpolation Equation:']);

Expanded Interpolation Equation:

disp(['P(x) = ' expanded_eqn]);

P(x) = 28.6763 + -17.4761 * (2 * cos(1* acos(x)) - 1)^1 + 3.8073 * (2 * cos(2* acos(x)) - 1)^2 + -0.17701 * (2 * cos(3* acos(x)) - 1)^3 + -0.17183 * (2 * cos(4* acos(x)) - 1)^4 + -1.7569 * (2 * cos(5* acos(x)) - 1)^5

% Plot the data
figure;
plot(x_data, y_data, 'ro', 'MarkerSize', 10); % plot data points
hold on;
plot(x_query*(x_max-x_min)/2+(x_max+x_min)/2, y_interp, 'b-', 'LineWidth', 2); % plot interpolating polynomial
xlabel('Year');
ylabel('CO2 Concentration');
title('Chebyshev Interpolation');
legend('Given Data', 'Interpolating Polynomial', 'Location', 'best');
grid on;

%% Estimate P(x) for specific x values
% Evaluate the interpolating polynomial at specific x-values
x_values = [1915, 2010, 2050];
y_values = zeros(size(x_values));

% Evaluate the interpolating polynomial at the given x-values
for i = 1:length(x_values)
    x_norm = 2 * (x_values(i) - x_min) / (x_max - x_min) - 1;
    for j = 1:n
        y_values(i) = y_values(i) + a(j) * T(j, x_norm);
    end
end

disp(['For x = 1915, y = ', num2str(y_values(1))]);

For x = 1915, y = 27.4652

disp(['For x = 2010, y = ', num2str(y_values(2))]);

For x = 2010, y = 16.5339

disp(['For x = 2050, y = ', num2str(y_values(3))]);

For x = 2050, y = -30.148

toc(t)

Elapsed time is 0.692066 seconds.

Answer 1

[Mathieu NOE]

hello

maybe this ?

%% Data input
x_data = [1830; 1920; 1940; 1975; 2000; 2020]; %x values
y_data = [0.089; 3.52; 4.86; 17.05; 25.5; 35.01]; %y values

% Chebyshev  polynomial interpolation
    x = linspace(x_data(1),x_data(end),25);

    % Barycentric interpolation
    m = numel(x_data);
    ts = -1+1/m:2/m:1-1/m;
    c = (-1).^(0:m-1).*sin(pi*(ts+1)/2);
    numer = zeros(size(x));
    denom = zeros(size(x));
    exact = zeros(size(x));

    for j = 1:m
        xdiff = x-x_data(j);
        temp = c(j)./xdiff;
        numer = numer + temp*y_data(j);
        denom = denom + temp;
        exact(xdiff==0) = j;
    end

    y = numer./denom;
    jj = find(exact);
    y(jj) = y_data(exact(jj));

% plot
plot(x_data,y_data,'*b',x,y,'r');

[Mathieu NOE]

hello again

problem solved ?

[Ignatius Tommy Pratama]

Hello Sir,

Thank you very much for your help. The problem is solved.

Again, thank you for your help. Have a great day.

[Mathieu NOE]

Tx

as always, my pleasure !