diff --git a/Function/Orbit_Propagation.m b/Function/Orbit_Propagation.m
new file mode 100644
index 0000000000000000000000000000000000000000..3ee3b9dfdbfbfe255f8da2f38ca55b3f20b8439f
--- /dev/null
+++ b/Function/Orbit_Propagation.m
@@ -0,0 +1,62 @@
+%% Function Description (Orbit Propagation)
+% The function `Orbit_Propagation` calculates the time derivatives of the state vector for a spacecraft
+% or celestial body undergoing two-body orbital dynamics. It is designed to be used with MATLAB's numerical
+% integrators like `ode45` to simulate orbital trajectories.
+%
+% This function takes in the current time `t`, the current state vector `y`, and the gravitational parameter
+% `mu` of the central body. The state vector `y` contains the position vector `r` and the velocity vector `v`,
+% both of which are 3-dimensional. The function computes the norm of the position vector `r_norm` to calculate
+% the acceleration due to gravity, which is directed towards the central body and inversely proportional to the
+% square of the distance from the central body.
+%
+% The differential equations that describe the motion are:
+% - `drdt`: The rate of change of the position, which is simply the velocity `v`.
+% - `dvdt`: The rate of change of the velocity, which represents the gravitational acceleration and is calculated as
+%   `(-mu * r) / r_norm^3`. This term ensures that the acceleration vector is always pointing towards the central
+%   body and is decreasing with the cube of the distance, as per Newton's law of universal gravitation.
+%
+% The output of the function is `dydt`, a vector that concatenates `drdt` and `dvdt`, thus providing the derivatives
+% of both position and velocity that are needed for the integration process.
+%
+% Usage:
+% dydt = Orbit_Propagation(t, y, mu)
+%
+%
+% Inputs:
+% t - Current time (not used in the calculation as the two-body problem is autonomous)
+% y - Current state vector, [position; velocity]
+% mu - Gravitational parameter of the central body (km^3/s^2)
+%
+% Output:
+% dydt - Derivative of the state vector, [velocity; acceleration]
+%
+% Author: Thomas West
+% Date: April 18, 2024
+
+%% Function
+% Define the differential equations function for the two-body problem
+function dydt = Orbit_Propagation(t, y, mu)
+    % Unpack the current state
+    % 'r' extracts the first three elements from the state vector 'y', representing the position vector.
+    % 'v' extracts elements four to six from the state vector 'y', representing the velocity vector.
+    r = y(1:3);
+    v = y(4:6);
+    
+    % Norm of the position vector
+    % 'r_norm' calculates the norm (magnitude) of the position vector 'r'.
+    % This is used in the denominator for the gravitational force calculation to normalize the position vector.
+    r_norm = norm(r);
+    
+    % Equations of motion for the two-body problem
+    % 'drdt' represents the rate of change of position, which is equal to the velocity vector 'v'.
+    drdt = v;
+    % 'dvdt' represents the rate of change of velocity, which is the gravitational acceleration.
+    % This acceleration is calculated as `(-mu * r) / r_norm^3`, where 'mu' is the gravitational parameter of the central body.
+    % The term 'r_norm^3' scales the gravitational force by the cube of the distance, according to Newton's law of universal gravitation.
+    dvdt = -mu * r / r_norm^3;
+    
+    % Output the derivative of the state
+    % 'dydt' is the output vector that concatenates 'drdt' and 'dvdt'.
+    % This vector provides the derivatives of both position and velocity, which are needed for numerical integration.
+    dydt = [drdt; dvdt];
+end