function qp_filter
omega = 10; beta = 20;
L = 1;
Omega = 1.1;
initial_w = [0.005,0.005,0,0,0,0.01];
[times,sols] = ode45(@(t,w) quadrupole(t,w,Omega,omega,beta),[0:1:120],initial_w);
animate_quadrupole(times,sols)

end
    function dwdt = quadrupole(t,w,Omega,omega,beta)
       x=w(1); y=w(2); z=w(3); vx=w(4); vy=w(5); vz=w(6);
       dxdt = vx; dydt=vy; dzdt=vz;
       dvxdt=Omega^2*(1+beta*cos(omega*t))*x;
       dvydt=Omega^2*(1+beta*cos(omega*t))*y;
       dvzdt=0;
       dwdt=[dxdt;dydt;dzdt;dvxdt;dvydt;dvzdt];
    end
    function [eval,stop_flag,edir] = event_function(t,w,L)
        x=w(1); y=w(2); z=w(3); vx=w(4); vy=w(5); vz=w(6);
        eval = L-z;
        stop_flag = 1;
        edir = -1;
    end

function animate_quadrupole(times,sols)

    [XX,YY,ZZ] = cylinder(0.1);
    XX1 = XX - 0.15;
    YY1 = YY - 0.15;
    XX2 = XX + 0.15;
    YY2 = YY + 0.15;
    figure1 = figure;
    for i=1:length(times)
        clf
     axes1 = axes('Parent',figure1,'DataAspectRatio',[1 1 1]);

     view(axes1,[-79.5 22]);
     grid(axes1,'on');
     xlim([-0.6 0.6]);
     ylim([-0.6 0.6]);
     zlim([0 1.2]);
     hold on
     surf(XX1,YY1,ZZ)
     surf(XX2,YY2,ZZ)
     surf(XX1,YY2,ZZ)
     surf(XX2,YY1,ZZ)      
     plot3(sols(1:i,1),sols(1:i,2),sols(1:i,3),...
         'LineWidth',2,'Color',[0 1 0]);
     plot3(sols(i,1),sols(i,2),sols(i,3),...
         'ro','MarkerSize',11,'MarkerFaceColor','r');
     pause(0.001)
    end
end
