//Феномен Рунге
function y=my_func(x)
y = 1./(1 + x.^2)
endfunction

function y=newton_method(xCount, x, xInterpolCount, Xn, Yn)
//заполняем нулями и высчитываем диагональ
    for i=1:xInterpolCount
        for j=1:xInterpolCount
            if i==j then
                M(i,j) = Yn(i);
            else
                M(i,j) = 0;
            end;
        end;
    end;

//высчитываем разделённые разности
    for i=1:+1:xInterpolCount-1
        for j=1:xInterpolCount-i
            M(j,j+i)=(M(j+1,j+i) - M(j,j+i-1)) / (Xn(j+i)-Xn(j));
        end
    end;

    for i = 1:xCount
        y(i) = Yn(1);
        if (xInterpolCount > 1)
        then
        p=1;
            for u=2:+1:xInterpolCount
                p = p*(x(i) - Xn(u-1));
                y(i) = y(i) + (M(1,u)*p);
            end
        end
    end
endfunction

function p=PolyLagr(m,n,x,xt,ft)
    // m кол-во Х  n- кол-во узлов
    //x-наши задаваемые Х     xt-узлы
    // ft - знач в узлах
    for k=1:m
        sum=0
        for i=1:n
            w(i)=1
            for j=1:n
                if (j ~= i)
            w (i) =w(i) * (x (k) - xt (j)) / (xt (i) - xt (j));
                end
            end
        sum = sum + ft(i) * w(i);
    end
    p(k) =sum;
    end
endfunction

function RP=RealP(n,f,f1)
    // f-данная функция   f1- интерполяц. полином
    for i=1:n
    RP(i)= abs(f(i)-f1(i))
    end
endfunction

// M- максимум в 10 производной
function Pr=Proizv(x,n) //n - степень

   i=0
   for  j=1:2:n //четные
       d1(j)=%i*((-1)^(j+i))*(factorial (j-1)./(1+%i*x)^(j+1)
       d2(j)=%i*((-1)^(j+i+1))*(factorial (j-1)./(1-%i*x)^(j+1)
       i=i+1
   end
   c=0
   for
       j=2:2:n //нечетное
       d1(j)=((-1)^(j+c))*(factorial (j-1)./(1+%i*x)^(j+1)
       d2(j)=((-1)^(j+c+1))*(factorial (j-1)./(1-%i*x)^(j+1)

     c=c+1
   end
   for i=1:n
   Pr(i)=0.5*d1(i)+0.5*d2(i)
   end
endfunction

function ()
    i=1:
    endfunction


xCount = 100;
xInterpolCount = 10;
a = -5;
b = 5;
h1 = (b-a)/(xCount - 1); //шаг для иксов функции
h2 = (b-a)/(xInterpolCount - 1); //шаг для узлов функции

x1 = [a:h1:b]; //иксы функции
x2 = [a:h2:b]; //узлы интерполяции

f1 = my_func(x1); //функция в иксах
F1 = my_func(x2); //функция в узловых точках
fN = newton_method(xCount, x1, xInterpolCount, x2, F1);//полином Ньютона
fL=PolyLagr(xCount,xInterpolCount,x1,x2,F1)
RealN=RealP(xCount,f1,fN)
RealL=RealP(xCount,f1,fL)
subplot(3,4,1)
plot(x1,f1,"r",x1,fN,"b");
xtitle('кр.-исходная функция, Син.-Ньютон')

subplot(3,4,2)
plot(x1,f1,"r",x1,fL,"g");
xtitle('кр.-исходная функция, зел.-Лагранж')

subplot(3,4,3)
plot(x1,RealN,"g");
xtitle('реал. погр. пол.Ньютона')

subplot(3,4,4)
plot(x1,RealL,'r');
xtitle('реал. погр. пол.Лагранжа')