Is 2000 materials a lot? Is there a better way?
Is 2000 materials a lot? It's a big number, but Blender can handle it. The real question is; is there a more manageable way. Maybe a way that lets you define one material, and use some property of the individual objects as a source for the color?

Yes there is. There's two parts to this answer.

Step 1 : Vertex Colors (layer)

Using a obj.data.vertex_colors on the Object's Mesh to set the RGB of the faces. This is per vertex per face, the code is very simple.

import bpy

obj = bpy.data.objects['Cube']
mesh = obj.data
polygons = mesh.polygons
vcols = mesh.vertex_colors

if len(vcols) == 0:
    vcols.new()

color_layer = vcols["Col"]

r, g, b = 0.3, 0.2, 0.9

for face in polygons:
    for v_idx in polygons[face.index].loop_indices:
        color_layer.data[v_idx].color = r, g, b




Step 2 : Shader Tree 'Attribute Node'

This Node has a text field where you can specify the vertex_colors layer by name, by default it's 'Col', but you can name it what you like.



To see the effect of the vertex_color layer in the 3d viewport (without rendering) Swith to Vertex Paint mode





Assigning one material to mutiple Objects

Once you have a material configured and want to share it over all objects in question, you can iterate over all these objects and do:

common_material = bpy.data.materials[some_named_material]
for obj in objects_in_question:
    obj.active_material = common_material




Things to consider

You raised another question about whether it might be easier to create the cubes as part of one object. That's totally doable if you keep track of the vertices and indices of the individual cubes ( at that point they are disjoint submeshes ) and when you know the face indices of each cube setting the vertex_color layer information is easy enough.


Doing it all as one object means you will need only 1 vertex color layer, and 1 material.
You could even use the foreach_set construct, it's perfect for meshes with uniform vertex-per-face counts (ie, all Tris, all Quads.. all N-gons)  

 # flattened_list = [r,g,b,r,g,b,r,g,b,r,g,b,.....]
 obj.data.vertex_colors[0].data.foreach_set("color", flattened_list)  


https://gist.github.com/zeffii/1a7d6800559f9ee902f8 shows how to do it for an Icosphere

import bpy

# default 80 faces, all Tris
obj = bpy.data.objects['Icosphere']  

# note the range is 80 * 3
flattened_list = []
[flattened_list.extend([0.5, 0.3, 0.7]) for i in range(80*3)]

obj.data.vertex_colors[0].data.foreach_set("color", flattened_list) 
obj.data.update()

Keeping the objects separate means you must use a vertex_color layer for each object, but you can share the common material.
Using thousands of materials for thousands of objects is certainly OK too.


Which of these options is the best depends on your scenario.


Code simplicity vs Efficiency
Efficiency vs 'How important is Efficiency'


How important is Efficiency: Is a High real-time framerate a requirement or is it being rendered and is the efficiency of your method not a real concern.