f = @(x) 3.*cos(1.5.*x).*exp(x/4)+4.*sin(3.5.*x).*exp(-3.*x)+4.*x;
a = 2.5;
b = 3.3;
alpha = 2/3;
betta = 0;

p= @(x)((x - a).^(-alpha)).*((b - x).^(-betta));
x = [];
mu = [];
for i = 0:1:2
    mu = [mu; integral(@(x)p(x).*x.^i,a,b)];
    x = [x;linspace(a,b,3).^i];
end
A = x\mu;
q = integral(matlabFunction(sym(f).*sym(p)),a,b);
q = abs(q - f(linspace(a,b,3))*A);
func_der = matlabFunction((-1).*abs(diff(sym(f), 3)));
M = (-1).* max([func_der(fminunc(func_der, a)), func_der(fminunc(func_der, (a+b)/2)), func_der(fminunc(func_der, b))]);
X = linspace(a,b,3);
omega = @(x) (x-X(1)).*(x-X(2)).*(x-X(3));
R = (M/factorial(3))*integral(matlabFunction(abs(sym(p)*sym(omega))), a, b)