Revision history [back]

In the end, I used ode_solver and an approximation of the delta function to get things working. Here is the code:

@interact(layout=[['A','p','interval'],['alpha','beta','omega']])
def lab(A=input_box(2,label='A',width=15),p=input_box(2,label='p',width=15),alpha = input_box(1,label='alpha',width=15),beta = input_box(1/sqrt(3),label='beta',width=15),omega = input_box(1,label='omega',width=15),interval=input_box((0,6),label='interval',width=20) ):
    def f(q):
        w=0.001
        if q>p-w/2 and q<p+w/2:
            return(1/w)
        else:
            return 0

    def dy(t,y):
        return [ 1, y[2], -omega^2*y[1]+A*f(y[0]) ]

    T=ode_solver()
    T.function=dy
    T.ode_solve(y_0=[0,alpha,beta],t_span=interval,num_points=500)
    ans=T.solution
    ans=[a[1] for a in ans]
    ty=[[a[0],a[1]] for a in ans]
    tv=[[a[0],a[2]] for a in ans]
    p=line(ty,legend_label='y (position)')+line(tv,color='red',legend_label='v (velocity)')