READ,'Enter the initial circuit voltage (V): ',v0
READ,'Enter the final circuit voltage (V): ',vf
READ,'Enter the circuit capacitance (microfarads): ',c
c=c*1.0E-6
READ,'Enter the circuit resistance (ohms): ',r

ERASE

tau=r*c
dt=0.01*tau
ntau=7
n=100*ntau
t=FINDGEN(n)*dt

tstring=STRING(v0,vf,c,r,tau $
       ,FORMAT='("V!i0!n=",I0," V, V!if!n=",I0," V, C=",E8.1," F,' $
              +' R=",E8.1," !4X!x, !4s!x=",E8.1," s")')

q=FLTARR(n)
q(0)=c*v0
v=q
v(0)=v0
i=v
i(0)=(vf-v0)/r

FOR k=1,n-1 DO BEGIN
    q(k)=q(k-1)+i(k-1)*dt
    v(k)=q(k)/c
    i(k)=(vf-v(k))/r
    ENDFOR

i=ABS(i)

tmax=MAX(t)
qmax=MAX(q,MIN=qmin)
vmax=MAX(v,MIN=vmin)
imax=MAX(i,MIN=imin)

!X.MARGIN=[13,3]

!P.MULTI=[0,1,3,0,0]
multiplot,/init

multiplot
PLOT,t,q,XRANGE=[0,tmax],YRANGE=[qmin,qmax],YTITLE='Charge (C)',TITLE=tstring

multiplot
PLOT,t,v,XRANGE=[0,tmax],YRANGE=[vmin,vmax],YTITLE='Voltage (V)'

multiplot
PLOT,t,i,XRANGE=[0,tmax],YRANGE=[imin,imax],YTITLE='Current (A)',XTITLE='Time (s)'

multiplot,/reset
!P.MULTI=0

END