clear all;
close all;
clc;

opt = input('Optiune: 1 = primul set; 2 = al doilea set; ');
if opt == 1
    xi=[-1 -0.96 -0.86 -0.79 0.22 0.5 0.93];
    yi=[-1 -0.151 0.894 0.986 0.895 0.5 -0.306];
    heval = 0.05;
    ni = 7;
else
    xi = [3 5.4 8.8 12.2 14 17.8 20.6 21.3 23.6 26];
    yi = [1540 1535 1530 1525 1520 1515 1510 1505 1500 1495];
    ni = 10;
    heval = 0.5;
end

xe = (xi(1):heval:xi(ni))';
npev = length(xe);
if (xe(npev) < xi(ni))
    npev = npev+1;
    xe = [xe; xi(ni)];
end
npev
mg = input('Introdu gradul polinomului de interpolare 1<=mg<=ni-1: ');
ye = zeros(npev,1);
for ip = 1:npev
    ye(ip) = INTPL(ni,xi,yi,mg,xe(ip));
    fprintf('xe(%g)=%g   ye(%g)=%g; \n',ip,xe(ip),ip,ye(ip));
end
if opt == 1
    plot(xi, yi, 'or', xe, ye, '-b');
    grid on;
else
    p = polyfit(xi, yi, 1);
    ye = polyval(p, xe);
    plot(xi, yi, 'or', xe, ye, '-b');
    grid on;
end


function ye = INTPL(ni,xi,yi,mg,xe)
    i=1;
    while ((i<=ni) && (xe>xi(i)))
        i = i+1;
    end
    Indmax = fix(i+mg/2);
    if (Indmax <= mg)
        Indmax = mg+i;
    end
    if (Indmax > ni)
        Indmax = ni;
    end
    Indmin = Indmax-mg;
    a=1;
    for i = Indmin:Indmax
        if (xe == xi(i))
            ye = yi(i);
            return;
        end
        a = a*(xe-xi(i));
    end
    ye=0;
    for i = Indmin:Indmax
        PL = yi(i)*a/(xe-xi(i));
        for j = Indmin:Indmax
            if (i ~= j)
                PL = PL/(xi(i)-xi(j));
            end
        end
        ye = ye+PL;
    end
end