how to fit quadratic with plateau condition in matlab?

Question

The values of x and y are as follows:

x=[0,80,100,120,150], 
y=[1865,3608,4057,4343,4389]

My equations are as follows:

y=a(x)^2+b(x)+c, x < xc(Ymax) ....(1)

y=yp, x >= xc(Ymax) ....(2)

I want to fit these two equations, where first equation will provide a quadratic fit whereas the second equation will provide a linear one. The problem is when we are trying to fit, it provide us best fit in quadratic form (eq.1), it did not include the linear fit (eq.2).

We are also providing the MATLAB code that we have used:

yf = @(b,x) b(1).*x.^2+b(2)*x+b(3);
B0 = [0.006; 21; 1878];
[Bm,normresm] = fminsearch(@(b) norm(y - yf(b,x)), B0);
a=Bm(1);
b=Bm(2);
c=Bm(3);
xc=-b/(2*a);
p=c-(b^2/(4*a));
 if (x < xc) 
      yfit = a.*x.^2+ b*x+c;
 else yfit = p;
 end
plot(x,yfit)
hold on 
plot(x,y,'*')
hold off

We are providing the graph (fig(1)) which we got from our results. For further understanding, we are also proving a graph (fig(2)) which shows how the result should come.

Answer 1

I modified the input data so you can see the result:

%input data
x = [0;     20;     50;     80;     100;    120;    150;    160;    180;    200;    220;     240]; 
y = [1865;  2500;  3000;   3608;   4057;   4343;   4389;    4250;   4300;   4280;   4350;    4200];

%find the index of the maximum y value
[~, idx] = max(y);

%initial parameter values
a0 = -1;
b0 = 1;
c0 = 2;

p0 = [a0 b0 c0];

%optimization options
options = optimset('Display','iter');

%optimization call
p = fminsearch(@fun_so_180730, p0, options, x, y, idx);

%plot of the fitted function based on two parts (quadratic and linear)
x_fit_1 = (0:2:x(idx))'; 
x_fit_2 = (x(idx):2:max(x))';

y_fit_1 = p(1)*x_fit_1.^2 + p(2)*x_fit_1 + p(3);
y_fit_2 = ones(size(x_fit_2))*y_fit_1(end);

figure;
plot(x, y, 'bo');
hold on;
plot([x_fit_1; x_fit_2], [y_fit_1; y_fit_2], 'r-');
hold off;
grid on;

The optimization function as a separate file:

function [f] = fun_so_180730(p, x, y, idx)

    a = p(1);
    b = p(2);
    c = p(3);

    % calculation of the current function
    y_prime = zeros(size(y));

    %quadratic part
    y_prime(1:idx) = a*x(1:idx).^2 + b*x(1:idx) + c;

    %linear part if there are elements after the maximum y value
    if (idx < numel(y))
        y_prime(idx+1:end) = y_prime(idx);
    end

    % current cost as a vector norm
    f = norm(y - y_prime);
end

Result: