I am at a bit of a loss as to understanding what you are doing and what you want to do. If you want a smoother result from your ODE integration, my first approach would be to use a finer increment for ‘tspan’, so instead of:
tspan=[0:0.1:20];
use:
tspan=linspace(0, 20, 500);
I don’t quite follow the filter design, though. What do you want to get from filtering your signal? Did you use the freqz function with it so you have some idea of its frequency and phase chatacteristics? There may be more appropriate filter designs for what you want to do.
