Organic yet accurate modeling with the golden spiral
How about a fully math based approach, here's an interesting paper on scipress.org called Growth in Plants: A Study in Numbers by Jay Kappraff.
from the paper:

Keywords: Phyllotaxis, Fibonacci Numbers, Golden Mean, Continued Fractions, Farey
Series, Wythoff’s Game, Zeckendorf Notation

Blender Addon(s)
There is an addon called Add 3d Function Surface which can be coerced into producing many interesting shapes. It is included in the Extra Objects addon as &quot;Math Functions&quot;

If the purely function based approach fails you then you still have a variety of options. Some of which are a bit more interactive.
Modifiers and Scripting
You could model one part of the object and use a script to rotate it around and scale it up/down according to whatever algorithm you can imagine. For example a 2d version of a plant  on geometry daily, it's simplicity is deceptive. I tend to use live coding tools like tributary.io to figure out the math using .js with immediate feedback. Usually if I can figure something out in 2D, the move to Blender and adding a third dimension isn't that much extra work.
An Array Modifier and an empty can give very similar results to the broccoli.

This is a cone, tilted slightly and shifted away from its origin, then added an Array modifier and used a tilted (rotated) and slightly scaled empty as an array object.
The fractal nature of the brocolli calls for daisy-chained array modifiers, each building on the complexity of the previous.

Elaboration
perhaps a short walk-through script might better convey what otherwise would become a wordy explanation. This script is merely intended to show what daisy-chained Array Modifiers can do, and shouldn't be perceived as nice python or nice code. The script should result in something like the sunflower. Run the code in an empty blend file with the 3D cursor at the centre of the scene.
import bpy
import math
import mathutils
from mathutils import Euler
# make sure 3d cursor is at the center of the scene

num_items_in_spiral = 13
num_spiral_arms = 29
z_rotation = 2.0 * math.pi / num_spiral_arms

# add cone, modify the mesh,
bpy.ops.mesh.primitive_cone_add(vertices=12, radius1=0.354406, depth=0.901211)
obj = bpy.context.active_object
obj.data.transform(mathutils.Matrix.Translation((-0.5, 0, 0)))

# add empty for spiral arm
bpy.ops.object.empty_add()
obj = bpy.context.active_object
obj.name = &quot;Empty_spiral&quot;
obj.location = (-0.8, 0, 0)
obj.rotation_euler = Euler((0, 0, 0.179609), 'XYZ')

# add empty for radial array
bpy.ops.object.empty_add()
obj = bpy.context.active_object
obj.name = &quot;Empty_radial_array&quot;
obj.rotation_euler = Euler((0, 0, z_rotation), 'XYZ')

# add modifiers to the cone
obj = bpy.data.objects[&quot;Cone&quot;]

mod = obj.modifiers.new(name=&quot;Array_1&quot;, type='ARRAY')
mod.use_relative_offset = False
mod.use_object_offset = True
mod.offset_object = bpy.data.objects[&quot;Empty_spiral&quot;]
mod.count = num_items_in_spiral

mod = obj.modifiers.new(name=&quot;Array_2&quot;, type='ARRAY')
mod.use_relative_offset = False
mod.use_object_offset = True
mod.offset_object = bpy.data.objects[&quot;Empty_radial_array&quot;]
mod.count = num_spiral_arms


The red cone is the primary object, the curve/spiral coming off it is the result of one array modifier on a translated and rotated Empty. The rest of the spirals are an extra array modifier on a different Empty (rotated around the scene origin by 2*math.pi / num_spiral_arms)
Whatever object/picture you choose to replicate, you can figure a lot out by drawing on top of the pictures and counting, but mostly experimenting. I think it's more accurate to consider these nice geometric objects as 3D fractals rather than the 3D embodiment of a golden ratio.