tspan = [0 4];

figure(1)
hold on
y0 = -4:0.1:4;    % List of initial conditions
for i = 1:length(y0)
    [t,y] = ode45('ic_variation_rhs', tspan, y0(i));
    plot(t,y)
end