It would help to have ‘dirfield’ (whatever it is). The other problem is that with one differential equation, ‘y0’ has to have one value, not a vector.
Put the plot and the differential equation solution in a loop, and loop over the ‘y0’ values —
f = @(t,y) (3880 - 0.817*y + (731000*y.^7.88)/(116000.^7.88 + y.^7.88));
% dirfield(f,0:10:100,0:1000:10000);
y0 = 0:100:8000;
figure
hold on
for k = 1:numel(y0)
[ts,ys] = ode45(f,[0,50],y0(k));
plot(ts,ys)
end
hold off
grid
xlabel('tx')
ylabel('ys')
I am not certain what you want as a result.
.
