You can use 'roots'. For each value of x you can compute the corresponding value of each of the four coefficients of the powers of y. For example the powers of y^3 can be collected to form a single coefficient which would be:
c3(x) = -12211.07+187137.42*x-935832.82*x^2+1384550.57*x^3
Use a for-loop:
x = linspace(0.15,0.39,n); % You choose n (= 25?)
y = zeros(n,3);
for k = 1:n
y(k,:) = roots([c3(x(k)),c2(x(k)),c1(x(k)),c0(x(k))]).'; % Get three y roots each trip
end
This gives you an n x 3 array of the corresponding possible values of y, some of which may be complex-valued, depending on the numbers you are using.
You may wish to reorder the rows in y you obtain this way in order to preserve continuity in the three branches. That would be an extra step.
