% Preset delta
delta = [0, 0.3, 0.8, 1.6, 5];
step = 0.01;

x = (0:(max(delta)/step)) * step;
f = log(1+exp(-x));
%p(x)=? 
p = ones(size(x));
pf = f .* p;

error = 1;
while error > 1e-6,
    %Calculate y
    for i = 1 : 4,
        first = floor(delta(i) / step) + 1;
        last = floor(delta(i+1) / step);
    
        y(i) = sum(pf(first:last)) / sum(p(first:last));
    end

    % recalc delta
    new_fd = (y(1:3) + y(2:4)) / 2;
    new_d = - log(exp(new_fd) - 1);

    error = sum((new_d - delta(2:4)).^2);
    delta(2:4) = new_d;
end
format long
y
new_d
error

x = 0:0.1:4;
plot(x, log(1+exp(-x)), 'r'); hold on;
xlabel('x');
ylabel('y');

xi = delta(1:4);
xp = [xi([1,2,2,3,3,4,4]), 4];
yp = y([1,1,2,2,3,3,4,4]);
plot(xp,yp, 'b--');

legend('f(x)', 'g(x)');