diff --git a/.gitignore b/.gitignore
index 3093fe03..81533fd7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -58,3 +58,6 @@ __pycache__
 # Custom config
 /custom.py
 /numdot_config.py
+
+# Program Debug Database
+*.pdb
\ No newline at end of file
diff --git a/demo/demos/boids_simulation/boid.png b/demo/demos/boids_simulation/boid.png
new file mode 100644
index 00000000..1a8d7c76
Binary files /dev/null and b/demo/demos/boids_simulation/boid.png differ
diff --git a/demo/demos/boids_simulation/boid.png.import b/demo/demos/boids_simulation/boid.png.import
new file mode 100644
index 00000000..63ca0840
--- /dev/null
+++ b/demo/demos/boids_simulation/boid.png.import
@@ -0,0 +1,34 @@
+[remap]
+
+importer="texture"
+type="CompressedTexture2D"
+uid="uid://wxlyngv1jx7x"
+path="res://.godot/imported/boid.png-3c7feec7ec9b604bb681557bcbd1345a.ctex"
+metadata={
+"vram_texture": false
+}
+
+[deps]
+
+source_file="res://demos/boids_simulation/boid.png"
+dest_files=["res://.godot/imported/boid.png-3c7feec7ec9b604bb681557bcbd1345a.ctex"]
+
+[params]
+
+compress/mode=0
+compress/high_quality=false
+compress/lossy_quality=0.7
+compress/hdr_compression=1
+compress/normal_map=0
+compress/channel_pack=0
+mipmaps/generate=false
+mipmaps/limit=-1
+roughness/mode=0
+roughness/src_normal=""
+process/fix_alpha_border=true
+process/premult_alpha=false
+process/normal_map_invert_y=false
+process/hdr_as_srgb=false
+process/hdr_clamp_exposure=false
+process/size_limit=0
+detect_3d/compress_to=1
diff --git a/demo/demos/boids_simulation/boid.tres b/demo/demos/boids_simulation/boid.tres
new file mode 100644
index 00000000..1cb8e4fb
--- /dev/null
+++ b/demo/demos/boids_simulation/boid.tres
@@ -0,0 +1,4 @@
+[gd_resource type="CompressedTexture2D" format=3 uid="uid://ci2unstco8670"]
+
+[resource]
+load_path = "res://.godot/imported/boid.png-3c7feec7ec9b604bb681557bcbd1345a.ctex"
diff --git a/demo/demos/boids_simulation/boid.tscn b/demo/demos/boids_simulation/boid.tscn
new file mode 100644
index 00000000..6edafb34
--- /dev/null
+++ b/demo/demos/boids_simulation/boid.tscn
@@ -0,0 +1,22 @@
+[gd_scene load_steps=2 format=3 uid="uid://b65qrclghjwb1"]
+
+[sub_resource type="CompressedTexture2D" id="CompressedTexture2D_nx014"]
+
+[node name="Boid" type="Control"]
+modulate = Color(0.1634, 0.76712, 0.86, 1)
+layout_mode = 3
+anchors_preset = 15
+anchor_right = 1.0
+anchor_bottom = 1.0
+grow_horizontal = 2
+grow_vertical = 2
+
+[node name="BoidTexture" type="TextureRect" parent="."]
+z_index = -1
+layout_mode = 0
+offset_left = -12.8
+offset_top = -25.6
+offset_right = 243.2
+offset_bottom = 358.4
+scale = Vector2(0.1, 0.1)
+texture = SubResource("CompressedTexture2D_nx014")
diff --git a/demo/demos/boids_simulation/boids_model.gd b/demo/demos/boids_simulation/boids_model.gd
new file mode 100644
index 00000000..b05b5608
--- /dev/null
+++ b/demo/demos/boids_simulation/boids_model.gd
@@ -0,0 +1,164 @@
+extends Node2D
+
+@export var boid_count: int = 20
+@export var speed: float = 200.0
+@export var range: float = 100.0
+@export var separation_weight: float = 1.0
+@export var alignment_weight: float = 0.5
+@export var cohesion_weight: float = 0.5
+
+@export var scale_factor: float = 0.1
+
+@export_category("Simulation parameters")
+@export var solver: BoidsSolver
+
+var texture = preload("res://demos/boids_simulation/boid.tres")
+
+var frame_time: float = 0.
+
+func _ready() -> void:
+	"""
+	Synchronize initial values with GUI, instantiate boid nodes, and initialize solver.
+	"""
+	# Synchronize initial values with GUI
+	boid_count = %NumberOfBoidsSlider.value
+	speed = %SpeedSlider.value
+	range = %RangeSlider.value
+	separation_weight = %SeparationSlider.value
+	alignment_weight = %AlignmentSlider.value
+	cohesion_weight = %CohesionSlider.value
+
+	# Instantiate Boid nodes under $Boids
+	if not has_node("Boids"):
+		var boids_container = Node2D.new()
+		boids_container.name = "Boids"
+		add_child(boids_container)
+	initialize_boids(boid_count)
+
+	solver.initialize()
+
+
+func _process(delta: float) -> void:
+	"""
+	Called every frame. Performs simulation step and updates GUI metrics.
+	
+	Parameters:
+	delta (float): Time since last frame in seconds.
+	"""
+	frame_time = Time.get_ticks_usec()
+	solver.simulation_step(delta)
+	frame_time = Time.get_ticks_usec() - frame_time
+
+	%FPSLabel.text = "FPS: " + str(Engine.get_frames_per_second())
+	%FrameTimeLabel.text = "Delta (ms): " + str(snappedf(frame_time/1000, 1e-3))
+
+
+func initialize_boids(target_count: int) -> void:
+	"""
+	Initializes boid nodes to match the specified target count.
+	
+	Parameters:
+	target_count (int): Desired number of boids.
+	"""
+	var current_count = $Boids.get_child_count()
+	if current_count > target_count:
+		var boids = $Boids.get_children()
+		for boid in range(target_count, current_count).map(func(index): return boids[index]):
+			boid.queue_free()
+	else:
+		for i in range(target_count-current_count):
+			add_boid(current_count+i)
+
+
+func add_boid(i: int):
+	"""
+	Adds a single boid sprite to the simulation.
+	
+	Parameters:
+	i (int): Index identifier for naming the new boid.
+	"""
+	var boid = Sprite2D.new()
+	boid.texture = texture
+	boid.set_modulate(Color.DARK_SLATE_GRAY)
+	boid.z_index = -1
+	boid.scale *= scale_factor
+	boid.name = "Boid" + str(i)
+	$Boids.add_child(boid)
+
+
+func _on_speed_slider_value_changed(value) -> void:
+	"""
+	Updates speed based on slider input and updates GUI label.
+	
+	Parameters:
+	value (float): New speed value from the slider.
+	"""
+	speed = value
+	%SpeedLabel.text = "Speed: " + str(speed)
+
+func _on_range_slider_value_changed(value) -> void:
+	"""
+	Updates range based on slider input and updates GUI label.
+	
+	Parameters:
+	value (float): New range value from the slider.
+	"""
+	range = value
+	%RangeLabel.text = "Range: " + str(range)
+
+func _on_separation_slider_value_changed(value) -> void:
+	"""
+	Updates separation weight based on slider input and updates GUI label.
+	
+	Parameters:
+	value (float): New separation weight value.
+	"""
+	separation_weight = value
+	%SeparationLabel.text = "Separation: " + str(separation_weight)
+
+func _on_alignment_slider_value_changed(value) -> void:
+	"""
+	Updates alignment weight based on slider input and updates GUI label.
+	
+	Parameters:
+	value (float): New alignment weight value.
+	"""
+	alignment_weight = value
+	%AlignmentLabel.text = "Alignment: " + str(alignment_weight)
+
+func _on_cohesion_slider_value_changed(value) -> void:
+	"""
+	Updates cohesion weight based on slider input and updates GUI label.
+	
+	Parameters:
+	value (float): New cohesion weight value.
+	"""
+	cohesion_weight = value
+	%CohesionLabel.text = "Cohesion: " + str(cohesion_weight)
+
+func _on_restart_button_pressed() -> void:
+	"""
+	Resets and reinitializes the solver upon pressing the restart button.
+	"""
+	solver.initialize()
+
+func _on_solver_option_item_selected(index: int) -> void:
+	"""
+	Switches solver based on selection from options.
+	
+	Parameters:
+	index (int): Index of the newly selected solver.
+	"""
+	solver = $Solvers.get_child(index)
+	solver.initialize()
+
+func _on_number_of_boids_slider_value_changed(value) -> void:
+	"""
+	Adjusts the number of boids based on slider input and updates GUI label.
+	
+	Parameters:
+	value (int): New desired count of boids.
+	"""
+	boid_count = value
+	%NumberOfBoids.text = "Boids: " + str(boid_count)
+	initialize_boids(boid_count)
diff --git a/demo/demos/boids_simulation/gd_solver.gd b/demo/demos/boids_simulation/gd_solver.gd
new file mode 100644
index 00000000..8272eeb0
--- /dev/null
+++ b/demo/demos/boids_simulation/gd_solver.gd
@@ -0,0 +1,134 @@
+extends BoidsSolver
+
+var positions: Array[Vector2]
+var directions: Array[Vector2]
+
+var screen_size: Vector2
+
+func initialize() -> void:
+	"""
+	Initializes the solver by setting the screen size
+	and generating initial position and direction arrays.
+	"""
+	screen_size = params.get_viewport_rect().size
+
+	# Initialize position and direction vector
+	positions = initialize_position_array(params.boid_count)
+	directions = initialize_direction_array(params.boid_count)
+
+
+# Helper function to create position direction vector with length
+func initialize_position_array(length: int) -> Array[Vector2]:
+	"""
+	Creates array of random positions within screen bounds.
+
+	Parameters:
+	length (int): Number of positions.
+
+	Returns:
+	Array[Vector2]: Random position vectors.
+	"""
+	# Initialize position Array with |length| Vector2s
+	# Values are random 2D-positions on screen
+	var positions_xy: Array[Vector2] = []
+	for i in range(length):
+		positions_xy.append(Vector2(randf()*screen_size.x, randf()*screen_size.y))
+	return positions_xy
+
+
+# Helper function to create random direction vector with length
+func initialize_direction_array(length: int) -> Array[Vector2]:
+	"""
+	Creates array of normalized random direction vectors.
+
+	Parameters:
+	length (int): Number of directions.
+
+	Returns:
+	Array[Vector2]: Random direction vectors.
+	"""
+	# Initialize direction Array with |length| Vector2s
+	# Values are normalized 2D-vectors with random angle
+	var directions_xy: Array[Vector2] = []
+	for i in range(length):
+		var angle = 2*PI*randf()
+		var direction = Vector2(cos(angle), sin(angle))
+		directions_xy.append(direction)
+	return directions_xy
+
+
+func simulation_step(delta: float) -> void:
+	"""
+	Updates positions and directions based on separation,
+	alignment, and cohesion during simulation.
+
+	Parameters:
+	delta (float): Time since last frame.
+	"""
+	# Check if boid_count has been changed, update vector sizes accordingly
+	var boid_count_difference = params.boid_count-positions.size()
+	if boid_count_difference < 0:
+		positions.resize(params.boid_count)
+		directions.resize(params.boid_count)
+	elif boid_count_difference > 0:
+		var new_positions := initialize_position_array(boid_count_difference)
+		var new_directions := initialize_direction_array(boid_count_difference)
+		positions.append_array(new_positions)
+		directions.append_array(new_directions)
+
+	# Calculate separation, alignment and cohesion per boid
+	var separations: Array[Vector2] = []
+	var alignments: Array[Vector2] = []
+	var cohesions: Array[Vector2] = []
+	for i in range(params.boid_count):
+		var separation := Vector2(0, 0)
+		var alignment := Vector2(0, 0)
+		var cohesion := Vector2(0, 0)
+		for j in range(params.boid_count):
+			var distance = (positions[i] - positions[j]).length()
+			if distance < params.range * 0.5 and distance != 0:
+				separation += (positions[i]-positions[j])/(distance**2)
+			if distance < params.range:
+				alignment += directions[j]
+				cohesion += positions[j]-positions[i]
+		if separation != Vector2(0, 0):
+			separation /= separation.length()
+		if cohesion != Vector2(0, 0):
+			cohesion /= cohesion.length()
+		if alignment != Vector2(0, 0):
+			alignment /= alignment.length()
+		separations.append(separation)
+		cohesions.append(cohesion)
+		alignments.append(alignment)
+
+	for i in range(params.boid_count):
+		# Apply separation, alignment and cohesion to direction vector and normalize
+		directions[i] += separations[i] * params.separation_weight * delta * 2.0
+		directions[i] += cohesions[i] * params.cohesion_weight * delta
+		directions[i] += alignments[i] * params.alignment_weight * delta
+		directions[i] /= directions[i].length()
+
+		# Move positions in directions by delta*speed
+		positions[i] += directions[i]*delta*params.speed
+
+		# Make boid positions wrap around at borders of screen
+		positions[i] = Vector2(fposmod(positions[i].x, screen_size.x), fposmod(positions[i].y, screen_size.y))
+	update_boids()
+
+
+func update_boids() -> void:
+	"""
+	Updates graphical positions and orientations of boids.
+	"""
+	var boids := params.get_node("Boids").get_children()
+	for i in range(params.boid_count):
+		var boid: Node2D = boids[i]
+
+		# Set position of boids by updating origin of transform
+		boid.transform.origin = positions[i]
+
+		# Set rotation of boids by aligning direction with up-vector of transform
+		var up := directions[i]
+		var right := Vector2(up.y, -up.x)
+		boid.transform.x = right*params.scale_factor
+		boid.transform.y = -up*params.scale_factor
diff --git a/demo/demos/boids_simulation/main.tscn b/demo/demos/boids_simulation/main.tscn
new file mode 100644
index 00000000..4c61bbf2
--- /dev/null
+++ b/demo/demos/boids_simulation/main.tscn
@@ -0,0 +1,200 @@
+[gd_scene load_steps=5 format=3 uid="uid://bujidut18snxp"]
+
+[ext_resource type="Script" path="res://demos/boids_simulation/boids_model.gd" id="1_2fw6l"]
+[ext_resource type="Script" path="res://demos/boids_simulation/gd_solver.gd" id="2_yi4li"]
+[ext_resource type="Script" path="res://demos/boids_simulation/nd_solver.gd" id="3_av8qh"]
+
+[sub_resource type="ImageTexture" id="ImageTexture_03ndn"]
+
+[node name="BoidsModel" type="Node2D" node_paths=PackedStringArray("solver")]
+script = ExtResource("1_2fw6l")
+solver = NodePath("Solvers/GDSolver")
+
+[node name="Solvers" type="Node" parent="."]
+
+[node name="GDSolver" type="Node" parent="Solvers" node_paths=PackedStringArray("params")]
+script = ExtResource("2_yi4li")
+params = NodePath("../..")
+
+[node name="NDSolver" type="Node" parent="Solvers" node_paths=PackedStringArray("params")]
+script = ExtResource("3_av8qh")
+params = NodePath("../..")
+
+[node name="ColorRect" type="ColorRect" parent="."]
+show_behind_parent = true
+z_index = -10
+custom_minimum_size = Vector2(1152, 648)
+anchors_preset = 5
+anchor_left = 0.5
+anchor_right = 0.5
+offset_right = 1152.0
+offset_bottom = 648.0
+grow_horizontal = 2
+color = Color(0.0431373, 0.0745098, 0.101961, 1)
+
+[node name="Labels" type="VBoxContainer" parent="."]
+offset_left = 23.0
+offset_top = 18.0
+offset_right = 223.0
+offset_bottom = 118.0
+metadata/_edit_group_ = true
+
+[node name="FPSLabel" type="Label" parent="Labels"]
+unique_name_in_owner = true
+layout_mode = 2
+theme_override_font_sizes/font_size = 30
+text = "FPS: 60"
+
+[node name="FrameTimeLabel" type="Label" parent="Labels"]
+unique_name_in_owner = true
+layout_mode = 2
+theme_override_font_sizes/font_size = 30
+text = "Frame time: 1"
+
+[node name="SliderOptions" type="VBoxContainer" parent="."]
+custom_minimum_size = Vector2(250, 0)
+offset_left = 22.0
+offset_top = 432.0
+offset_right = 272.0
+offset_bottom = 626.0
+metadata/_edit_group_ = true
+
+[node name="SeparationLabel" type="Label" parent="SliderOptions"]
+unique_name_in_owner = true
+layout_mode = 2
+theme_override_font_sizes/font_size = 30
+text = "Separation: 0.2"
+
+[node name="SeparationSlider" type="HSlider" parent="SliderOptions"]
+unique_name_in_owner = true
+layout_mode = 2
+max_value = 1.0
+step = 0.01
+value = 0.2
+
+[node name="AlignmentLabel" type="Label" parent="SliderOptions"]
+unique_name_in_owner = true
+layout_mode = 2
+theme_override_font_sizes/font_size = 30
+text = "Alignment: 0.4
+"
+
+[node name="AlignmentSlider" type="HSlider" parent="SliderOptions"]
+unique_name_in_owner = true
+layout_mode = 2
+max_value = 1.0
+step = 0.01
+value = 0.4
+
+[node name="CohesionLabel" type="Label" parent="SliderOptions"]
+unique_name_in_owner = true
+layout_mode = 2
+theme_override_font_sizes/font_size = 30
+text = "Cohesion: 0.6
+"
+
+[node name="CohesionSlider" type="HSlider" parent="SliderOptions"]
+unique_name_in_owner = true
+layout_mode = 2
+max_value = 1.0
+step = 0.01
+value = 0.6
+
+[node name="RestartButton" type="Button" parent="."]
+offset_left = 865.0
+offset_top = 79.0
+offset_right = 1140.0
+offset_bottom = 129.0
+theme_override_font_sizes/font_size = 30
+text = "Restart
+"
+
+[node name="SolverOption" type="OptionButton" parent="."]
+offset_left = 865.0
+offset_top = 23.0
+offset_right = 1137.0
+offset_bottom = 73.0
+theme_override_font_sizes/font_size = 30
+selected = 0
+item_count = 2
+popup/item_0/text = "GDScript"
+popup/item_1/text = "NumDot"
+popup/item_1/id = 1
+
+[node name="CanvasLayer" type="CanvasLayer" parent="."]
+
+[node name="TextureRect" type="TextureRect" parent="CanvasLayer"]
+texture_filter = 1
+offset_right = 40.0
+offset_bottom = 40.0
+size_flags_horizontal = 4
+size_flags_vertical = 4
+texture = SubResource("ImageTexture_03ndn")
+
+[node name="LegendInfected" type="TextureRect" parent="."]
+offset_left = 288.0
+offset_top = 541.0
+offset_right = 308.0
+offset_bottom = 561.0
+
+[node name="LegendRecovered" type="TextureRect" parent="."]
+offset_left = 288.0
+offset_top = 607.0
+offset_right = 308.0
+offset_bottom = 627.0
+
+[node name="SolverOptions" type="VBoxContainer" parent="."]
+custom_minimum_size = Vector2(200, 0)
+offset_left = 902.0
+offset_top = 432.0
+offset_right = 1130.0
+offset_bottom = 626.0
+metadata/_edit_group_ = true
+
+[node name="SpeedLabel" type="Label" parent="SolverOptions"]
+unique_name_in_owner = true
+layout_mode = 2
+theme_override_font_sizes/font_size = 30
+text = "Speed: 200"
+
+[node name="SpeedSlider" type="HSlider" parent="SolverOptions"]
+unique_name_in_owner = true
+layout_mode = 2
+max_value = 1000.0
+value = 200.0
+
+[node name="NumberOfBoids" type="Label" parent="SolverOptions"]
+unique_name_in_owner = true
+layout_mode = 2
+theme_override_font_sizes/font_size = 30
+text = "Boids: 50"
+
+[node name="NumberOfBoidsSlider" type="HSlider" parent="SolverOptions"]
+unique_name_in_owner = true
+layout_mode = 2
+min_value = 1.0
+max_value = 1000.0
+step = 0.5
+value = 50.0
+exp_edit = true
+
+[node name="RangeLabel" type="Label" parent="SolverOptions"]
+unique_name_in_owner = true
+layout_mode = 2
+theme_override_font_sizes/font_size = 30
+text = "Range: 100"
+
+[node name="RangeSlider" type="HSlider" parent="SolverOptions"]
+unique_name_in_owner = true
+layout_mode = 2
+max_value = 1000.0
+value = 100.0
+
+[connection signal="value_changed" from="SliderOptions/SeparationSlider" to="." method="_on_separation_slider_value_changed"]
+[connection signal="value_changed" from="SliderOptions/AlignmentSlider" to="." method="_on_alignment_slider_value_changed"]
+[connection signal="value_changed" from="SliderOptions/CohesionSlider" to="." method="_on_cohesion_slider_value_changed"]
+[connection signal="pressed" from="RestartButton" to="." method="_on_restart_button_pressed"]
+[connection signal="item_selected" from="SolverOption" to="." method="_on_solver_option_item_selected"]
+[connection signal="value_changed" from="SolverOptions/SpeedSlider" to="." method="_on_speed_slider_value_changed"]
+[connection signal="value_changed" from="SolverOptions/NumberOfBoidsSlider" to="." method="_on_number_of_boids_slider_value_changed"]
+[connection signal="value_changed" from="SolverOptions/RangeSlider" to="." method="_on_range_slider_value_changed"]
diff --git a/demo/demos/boids_simulation/metadata.json b/demo/demos/boids_simulation/metadata.json
new file mode 100644
index 00000000..1803fcd7
--- /dev/null
+++ b/demo/demos/boids_simulation/metadata.json
@@ -0,0 +1,6 @@
+{
+	"name": "Boids Simulation",
+	"description": "A simulation of flocking behavior modeled on how birds and other animals move collectively.",
+	"link": "https://en.wikipedia.org/wiki/Boids",
+	"authors": "versjon, kro-ma"
+}
diff --git a/demo/demos/boids_simulation/nd_solver.gd b/demo/demos/boids_simulation/nd_solver.gd
new file mode 100644
index 00000000..671ad341
--- /dev/null
+++ b/demo/demos/boids_simulation/nd_solver.gd
@@ -0,0 +1,174 @@
+extends BoidsSolver
+
+var positions: NDArray
+var directions: NDArray
+
+var rng := nd.default_rng()
+
+var screen_size: Vector2
+
+func initialize() -> void:
+	"""
+	Initializes solver by setting the screen size
+	and generating initial position and direction arrays.
+	"""
+	screen_size = params.get_viewport_rect().size
+
+	# Initialize position and direction vector
+	positions = initialize_position_array(params.boid_count)
+	directions = initialize_direction_array(params.boid_count)
+
+
+# Helper function to create position direction vector with length
+func initialize_position_array(length: int) -> NDArray:
+	"""
+	Generates an NDArray of random positions within screen bounds.
+
+	Parameters:
+	length (int): Number of positions to generate.
+
+	Returns:
+	NDArray: Array containing random position vectors.
+	"""
+	# Initialize position vector of shape [length, 2]
+	# Values are random 2D positions on screen
+	var positions_xy = rng.random([length, 2])
+	positions_xy.assign_multiply(positions_xy, nd.array([screen_size.x, screen_size.y], 2))
+	return positions_xy
+
+
+# Helper function to create random direction vector with length
+func initialize_direction_array(length: int) -> NDArray:
+	"""
+	Creates an NDArray of normalized random direction vectors.
+
+	Parameters:
+	length (int): Number of directions to generate.
+
+	Returns:
+	NDArray: Array containing random normalized direction vectors.
+	"""
+	# Initialize angle vector of shape [length]
+	# Values are random angles in [0, 2*PI), used to create direction vector
+	var angles := rng.random([length])
+	angles.assign_multiply(angles, 2.0 * PI)
+
+	# Initialize direction vector of shape [length, 2]
+	# Values are normalized 2D direction vectors according to angles
+	var directions_x := nd.cos(angles)
+	var directions_y := nd.sin(angles)
+	return nd.stack([directions_x, directions_y], 1)
+
+
+func simulation_step(delta: float) -> void:
+	"""
+	Executes a simulation step updating positions and directions
+	using NumDot operations for efficiency.
+
+	Parameters:
+	delta (float): Time elapsed since the last frame.
+	"""
+	# Check if boid_count has been changed, update vector sizes accordingly
+	var boid_count_difference = params.boid_count-positions.shape()[0]
+	if boid_count_difference < 0:
+		positions = positions.get(nd.range(params.boid_count), nd.range(2))
+		directions = directions.get(nd.range(params.boid_count), nd.range(2))
+	elif boid_count_difference > 0:
+		var new_positions := initialize_position_array(boid_count_difference)
+		var new_directions := initialize_direction_array(boid_count_difference)
+		positions = nd.vstack([positions, new_positions])
+		directions = nd.vstack([directions, new_directions])
+
+	# Move positions in directions by delta*speed
+	var offset := nd.multiply(directions, delta*params.speed)
+	positions.assign_add(positions, offset)
+
+	# Make boid positions wrap around at borders of screen
+	for axis in [0, 1]:
+		var positions_axis := positions.get(nd.range(params.boid_count), axis)
+		var wrap_positive := nd.greater(positions_axis, screen_size[axis]).as_type(nd.Int16)
+		var wrap_negative := nd.less(positions_axis, 0).as_type(nd.Int16)
+		wrap_positive.assign_multiply(wrap_positive, -screen_size[axis])
+		wrap_negative.assign_multiply(wrap_negative, screen_size[axis])
+		positions_axis.assign_add(positions_axis, wrap_positive)
+		positions_axis.assign_add(positions_axis, wrap_negative)
+		positions.set(positions_axis, nd.range(params.boid_count), axis)
+
+	# Create pair-wise position differences of boids with shape [n, n, 2]
+	var position_differences := nd.subtract(positions.get(null, &"newaxis", null), positions.get(&"newaxis", null, null))
+	# Calculate distances from position differences with shape [n, n]
+	var position_distances := nd.norm(position_differences, 2, 2)
+	# Mark every pair of boids with distance smaller than range in separation mask with shape [n, n]
+	var vision_mask := nd.less_equal(position_distances, params.range).as_type(nd.Int16)
+	# Mark every pair of boids with distance smaller than 0.5*range in separation mask with shape [n, n]
+	var separation_mask := nd.less_equal(position_distances, params.range*0.5).as_type(nd.Int16)
+
+	# Separation
+	# Calculate separation direction normalization divisor with shape [n, n]
+	var separation_normalization := nd.square(position_distances)
+	separation_normalization.assign_add(separation_normalization, nd.eye(params.boid_count))
+	# Normalize separation directions inversely proportional to distances with shape [n, n, 1]
+	var separation_directions := nd.divide(position_differences, separation_normalization.get(null, null, &"newaxis"))
+	# Ignore boids not marked in separation mask
+	separation_directions.assign_multiply(position_differences, separation_mask.get(null, null, &"newaxis"))
+	# Calculate sum of separation directions per boid with shape [n, 2]
+	var separations := nd.sum(separation_directions, 0)
+	# Make seperation directions point away from boids in separation range
+	separations.assign_multiply(separations, -1)
+	# Calculate separation direction normalization divisor with shape [n]
+	var separations_normalization := nd.norm(separations, 2, 1)
+	separations_normalization.assign_add(separations_normalization, nd.equal(separations_normalization, 0.0).as_type(nd.Int16))
+	# Normalize separation directions for each boid
+	separations.assign_divide(separations, separations_normalization.get(null, &"newaxis"))
+
+	# Alignment
+	# Ignore boids not marked in vision mask
+	var alignment_directions := nd.multiply(directions.get(null, &"newaxis", null), vision_mask.get(null, null, &"newaxis"))
+	# Calculate sum of alignment directions per boid with shape [n, 2]
+	var alignments := nd.sum(alignment_directions, 0)
+	# Normalize alignment directions for each boid
+	alignments.assign_divide(alignments, nd.norm(alignments, 2, 1).get(null, &"newaxis"))
+
+	# Cohesion
+	# Ignore boids not marked in cohesion mask
+	var cohesion_positions := nd.multiply(positions.get(null, &"newaxis", null), vision_mask.get(null, null, &"newaxis"))
+	# Calculate sum of cohesion positions with shape [n, 1]
+	var cohesions := nd.sum(cohesion_positions, 0)
+	# Find cohesion centers by calculating averages
+	cohesions.assign_divide(cohesions, nd.sum(vision_mask, 0).get(null, &"newaxis"))
+	# Calculate cohesion directions by taking difference between boids and respective cohesion centers
+	cohesions.assign_subtract(cohesions, positions)
+	# Calculate cohesion direction normalization divisor (dist to cohesion center if existing, else 1)
+	var cohesions_normalization = nd.norm(cohesions, 2, 1)
+	cohesions_normalization.assign_add(cohesions_normalization, nd.equal(cohesions_normalization, 0.0).as_type(nd.Int16))
+	# Normalize cohesion directions
+	cohesions.assign_divide(cohesions, cohesions_normalization.get(null, &"newaxis"))
+
+	# Update directions vector according to separation, alignment and cohesion with respective weights
+	directions.assign_add(directions, separations.assign_multiply(separations, params.separation_weight*delta*2))
+	directions.assign_add(directions, alignments.assign_multiply(alignments, params.alignment_weight*delta))
+	directions.assign_add(directions, cohesions.assign_multiply(cohesions, params.cohesion_weight*delta))
+
+	# Normalize direction vection lengths to 1
+	directions.assign_divide(directions, nd.norm(directions, 2, 1).get(null, &"newaxis"))
+
+	update_boids()
+
+
+func update_boids() -> void:
+	"""
+	Updates the graphical representations of boids
+	to match the computed positions and directions.
+	"""
+	var boids := params.get_node("Boids").get_children()
+	for i in range(min(params.boid_count, boids.size())):
+		var boid: Node2D = boids[i]
+
+		# Set position of boids by updating origin of transform
+		boid.transform.origin = positions.get_vector2(i, nd.range(2))
+
+		# Set rotation of boids by aligning direction with up-vector of transform
+		var up := directions.get_vector2(i, nd.range(2))
+		var right := Vector2(up.y, -up.x)
+		boid.transform.x = right*params.scale_factor
+		boid.transform.y = -up*params.scale_factor
diff --git a/demo/demos/boids_simulation/solver.gd b/demo/demos/boids_simulation/solver.gd
new file mode 100644
index 00000000..2fd207c1
--- /dev/null
+++ b/demo/demos/boids_simulation/solver.gd
@@ -0,0 +1,13 @@
+extends Node
+class_name BoidsSolver
+
+@export var params: Node2D
+
+func initialize() -> void:
+	pass
+
+func simulation_step(delta: float) -> void:
+	pass
+
+func update_boids() -> void:
+	pass