One of the easiest ways I can think of is to use a continuous-time filter, with initial conditions set to zero. The only caveat is the attenuation in the magnitude of the sine wave and a bit of shift in the absolute phase. The relative phase difference is still 72 degrees. The magnitude attenuation can be corrected by multiplying it with a factor as done in the following code.
This code requires the control system toolbox. The same can be done with ode45.
t = linspace(0, 3, 1000);
freq = 1;
phases = 0:deg2rad(72):deg2rad(72)*5;
y = sin(2*pi*freq*t + phases.');
filter_freq = 3;
fil = tf(2*pi*filter_freq, [1 2*pi*filter_freq]);
y_fil = zeros(size(y));
for i=1:size(y,1)
y_fil(i, :) = lsim(fil, y(i,:), t);
end
mul_factor = sqrt((2*pi*freq)^2+(2*pi*filter_freq)^2)/(2*pi*filter_freq);
y_fil = mul_factor*y_fil;
plot(t, y_fil)
