Do calculation for array of values - Python

I am trying to calculate and plot the amplitude amp of the solution of a time dependent t differential equation of motion (see rhs6) as a function of wd for multiple values of the force coefficient f found in f_array.

So far I have got the code working for plotting amp against wd for a single value of f. The result is a resonance peak:

The code for plotting amp against wd for one value of f (which produced the above image) is given below.

from pylab import *
from scipy.integrate import odeint
from numpy import *
import math

#Parameters.
k = 2.0
m = 1.0
w0 = (k/m)**(1/2)
alpha = 0.2
l = alpha/(2*m)
f = 1.0
wd = w0 + 0.025
beta = 0.2
t_fixed = 2.0

#Arrays.
t = linspace(0.0, 400.0, 400.0)
wd_array = linspace(w0-1.0, w0+1, 400.0)
f_array = linspace(10.0, 100.0, 3.0)

#Time step.
init_x = 0.0
init_v = 0.0
dx = 15.0
dv = 0.0
init_cond = [init_x,init_v]
init_cond2 = [init_x + dx,init_v + dv]

def rhs6(c,t,wd):
    c0dot = c[1]
    c1dot = -2*l*c[1] - w0*w0*c[0] + (f/m)*cos((wd)*t)
    return [c0dot, c1dot]

amp_array=[]
for wd in wd_array:
    res = odeint(rhs6, init_cond, t, args=(wd,))
    amp = max(res[:,0]) 
    amp_array.append(amp)

plot(wd_array, amp_array)
xlabel('Driving frequency, wd')
ylabel('Ampltiude, amp')    
show()

Now I wish to find amp against wd for multiple values of f. I have tried to do this by basically making a for loop statement over the f_array. However, my approach does not work and I get the error:

setting an element with a sequence.

As it's good to show an attempt, below is mine.

from pylab import *
from scipy.integrate import odeint
from numpy import *
import math

#Parameters.
k = 2.0
m = 1.0
w0 = (k/m)**(1/2)
alpha = 0.2
l = alpha/(2*m)
f = 1.0
wd = w0 + 0.025
beta = 0.2
t_fixed = 2.0

#Arrays.
t = linspace(0.0, 200.0, 200.0)
wd_array = linspace(w0-1.0, w0+1, 200.0)
f_array = linspace(10.0, 200.0, 3.0)

#Time step.
init_x = 0.0
init_v = 0.0
dx = 15.0
dv = 0.0
init_cond = [init_x,init_v]
init_cond2 = [init_x + dx,init_v + dv]

def rhs6(c,t,wd,f):
    c0dot = c[1]
    c1dot = -2*l*c[1] - w0*w0*c[0] + (f/m)*cos((wd)*t)
    return [c0dot, c1dot]

full_array = zeros(len(f_array))    
for index,force in enumerate(f_array):
    amp_list = []
    for wd in wd_array:
        res = odeint(rhs6, init_cond, t, args=(wd,force))
        amp = max(res[:,0]) 
        amp_list.append(amp)
        print(res)
    amp_array = array(amp_list)
    full_array[index] = amp_array

for f in full_array:
    plot(wd, amp)

show()

Any ideas?