tic;
samples = 10000000;
xBounds = LinearBounds(0,20);
yBounds = LinearBounds(-10,150000);
bounds = Bounds(xBounds,yBounds);
pointsInside = 0;
realValue = realIntegral(bounds.x);

for i = 0:samples
    toTestPoint = bounds.getRandomPoint;
    if getIsInside(toTestPoint)
        pointsInside = pointsInside+1;
    end
end


format longG
integral = (pointsInside/samples)*bounds.area
realValue
error = abs(realValue-integral)
samples
toc



function out = mainFunction(x)
    out = x^4-4*x^3+18*x^2-12*x-69;
end

function out = indefiniteIntegral(x)
    out = (1/5)* x^5 - x^4 + 6* x^3 - 6*x^2 - 69*x;
end

function out = getIsInside(p)
    arguments
        p Point
    end
    y_0 = mainFunction(p.x);
    out = abs(y_0) > abs(p.y) && y_0 * p.y >= 0;
end

function out = realIntegral(bounds)
    arguments
        bounds LinearBounds
    end
    out = indefiniteIntegral(bounds.higher) - indefiniteIntegral(bounds.lower);
end