Driven LCR circuit

from pylab import*
from math import*
c=1
l=1
r=1.5
q0=0
i0=0
v0=1
w=1.5
t=0.0
h=0.01
def f1(t,i,q):
    return((-1/l)*((q/c)+r*i-v0*sin(w*t)))
def f2(t,i,q):
    return(i)
i=[i0]
q=[q0]
time=[t]

while t<200:
    k1=h*f1(t,i0,q0)
    l1=h*f2(t,i0,q0)
    k2=h*f1(t + h/2.0,i0+k1/2.0,q0+l1/2.0 )
    l2=h*f2(t + h/2.0,i0+k1/2.0,q0+l1/2.0 )
    k3=h*f1(t + h/2.0,i0+k2/2.0,q0+l2/2.0 )
    l3=h*f2(t + h/2.0,i0+k2/2.0,q0+l2/2.0 )
    k4=h*f1(t + h, i0+k3, q0+l3 )
    l4=h*f2(t + h, i0+k3, q0+l3 )
      
    i1=i0+(k1 + 2*k2+ 2*k3 + k4)/6.0
    q1=q0+(l1 + 2*l2+ 2*l3 + l4)/6.0
    i.append(i1)
    q.append(q1)
    t=t+h
    time.append(t)
    i0=i1
    q0=q1

figure(1)
title("Current in driven LCR circuit")
xlabel(" Time")
ylabel("Current")
plot(time,i)
grid(True)

figure(2)
title(" Charge in driven LCR Circuit")
xlabel("Time")
ylabel("Charge")
plot(time,q)
grid(True)

#figure(3)
#title("Harmonic oscillator motion")
#xlabel("x")
#ylabel("Velocity")
#plot(i,q)
#grid(True)

show()