close all;
clc;
clear;

F='1000./(x.^2-9*x+68)'
func = inline(F);
F='(-2000*x+9000)./((x.^2-9*x+68).^2)'
DerivativeFunc=inline(F);
F='(((2000.*(2*x-9).^2)./(x.^2-9*x+68))-2000)./((x.^2-9*x+68).^2)'
Derivative2Func = inline(F);

simple =[];
modSimple =[];
RungeKutta =[];
modRungeKutta =[];
secondDerivativeOnePoint=[];
secondDerivativeTwoPoints=[]

for i= 10:100

		x = linspace(2, 4, i);
		[d, ind] = min(abs(x - 3));
		h=x(2)-x(1);
		k=rand * 0.001 - 0.0001;

		f12 = func(x(ind - 2)) + k;
    f11 = func(x(ind - 1)) + k;
    f0 = func(x(ind)) + k;
    f1 = func(x(ind + 1)) + k;
    f2 = func(x(ind + 2)) + k;
    f3 = func(x(ind + 3)) + k;



    p1 = (f1 - f0)/h;
    del1 = abs(DerivativeFunc(x(ind)) - p1);
    alpha1 = -log(del1)/log(ind);
    simple= [simple alpha1];


    p2 = (f1 - f11)/(2*h);
    del2 = abs(DerivativeFunc(x(ind)) - p2);
    alpha2 = -log(del2)/log(ind);
    modSimple= [modSimple alpha2];


    p3 = (-3 * f0 + 4 * f1 - f2)/(2*h);
    del3 = abs(DerivativeFunc(x(ind)) - p2);
    alpha3 = -log(del3)/log(ind);
    RungeKutta= [RungeKutta alpha3];


    p4 = (-11 * f0 + 18 * f1 - 9 * f2 + 2 * f3)/(6*h);
    del4 = abs(DerivativeFunc(x(ind)) - p4);
    alpha4 = -log(del4)/log(ind);
    modRungeKutta= [modRungeKutta alpha4];


    p5 = (f1 - 2 * f0 + f11)/(h^2);
    del5 = abs(Derivative2Func(x(ind)) - p5);
    alpha5 = -log(del5)/log(ind);
    secondDerivativeOnePoint= [secondDerivativeOnePoint alpha5];

    p6 = (- f12 + 16 * f11 - 30 * f0 + 16 * f1 - f2)/(12 * h^2);
    del6 = abs(Derivative2Func(x(ind)) - p6);
    alpha6 = -log(del6)/log(ind);
    secondDerivativeTwoPoints = [secondDerivativeTwoPoints alpha6];
end

i=10:100;
plot (i, simple, 'c');
hold on
plot (i, modSimple, 'g');
plot (i, RungeKutta, 'k');
plot (i, modRungeKutta, 'm');
plot (i, secondDerivativeOnePoint, 'r');
plot (i, secondDerivativeTwoPoints, 'y');
legend('Simpal', 'modSimple','RungeKutta', 'modRungeKutta', 'secondDerivativeOnePoint', 'secondDerivativeTwoPoints');
grid on;