diff --git a/Function/lambert_solver.m b/Function/lambert_solver.m
new file mode 100644
index 0000000000000000000000000000000000000000..f102b396b20f051ba536f641c364d5956fee9c98
--- /dev/null
+++ b/Function/lambert_solver.m
@@ -0,0 +1,163 @@
+%% Function Description (Lambert Solver)
+% This function, `lambert_solver`, calculates the initial and final velocity vectors (v0 and vf)
+% required to transfer a spacecraft between two orbits defined by position vectors r0 and rf in a given time (Dt),
+% under a central body with a gravitational parameter (mu). It employs Lambert's problem solution to determine
+% these velocities, which is a fundamental problem in orbital mechanics for designing transfer orbits.
+%
+% The function starts by validating the time of flight against feasible minimum and maximum values to ensure
+% it falls within a realistic range. It then calculates the cosine of the change in true anomaly (delta nu)
+% and uses it to compute a term called 'A', crucial for the successive approximation used to solve Lambert's problem.
+%
+% A loop is employed to iterate through possible values of the auxiliary variable 'psi' until the time of flight
+% calculated from the current guess is close enough to the desired time of flight, Dt. This iterative method
+% converges upon the correct psi value that aligns the computed time of flight with the given Dt.
+%
+% The function uses Stumpff functions, essential for dealing with universal variables in orbital mechanics,
+% to compute values needed during iterations. Upon convergence, the function calculates the Lagrange f and g
+% coefficients to derive the velocity vectors at the start and end of the transfer.
+%
+% The function returns:
+% - v0: Initial velocity vector at r0
+% - vf: Final velocity vector at rf
+% - exitflag: A flag indicating if the function converged (1) or not (0)
+%
+% Additionally, it checks the orbit transfer direction (DM) to decide if the transfer should be performed
+% in a prograde (short way) or retrograde (long way) direction relative to the motion.
+%
+% Usage:
+% [v0, vf, exitflag] = lambert_solver(r0, rf, Dt, mu, DM)
+%
+% Inputs:
+% r0 - Initial position vector
+% rf - Final position vector
+% Dt - Time of flight (seconds)
+% mu - Gravitational parameter of the central body (km^3/s^2)
+% DM - Direction of motion. DM = 1; % Transfer mode setting:
+% -1 for a Type II trajectory, which takes the longer path, more than 180 degrees around the central body.
+% This path may be chosen for specific mission requirements like thermal management or avoiding solar conjunctions.
+
+% +1 for a Type I trajectory, representing the shorter path, less than 180 degrees around the central body.
+% This is generally more energy-efficient and quicker, suitable for missions prioritizing speed and fuel efficiency.
+%
+% Outputs:
+% v0 - Initial velocity vector
+% vf - Final velocity vector
+% exitflag - Convergence flag (1 for success, 0 for failure)
+%
+% Author: Thomas West
+% Date: April 18, 2024
+
+%% Lambert Solver Function
+function [v0, vf, exitflag] = lambert_solver(r0, rf, Dt, mu, DM)
+ % Constants
+ tol = 1e-6; % Tolerance for numerical calculations
+ max_iter = 1000; % Maximum number of iterations to attempt for convergence
+ iter_count = 0; % Initialize iteration counter
+
+ % Constants for minimum and maximum transfer time in seconds
+ min_transfer_time = 1 * 86400; % Minimum time of flight, set to 1 days
+ max_transfer_time = 1000 * 86400; % Maximum time of flight, set to 1,000 days
+
+ % Pre-calculation check for time of flight
+ % Check if the provided time of flight (Dt) is within the feasible range.
+ if Dt < min_transfer_time || Dt > max_transfer_time
+ error('Time of flight is outside the feasible range (1 to 1,000 days).');
+ end
+
+ % Calculate the magnitudes of initial and final position vectors
+ r0_mag = norm(r0); % Magnitude of the initial position vector
+ rf_mag = norm(rf); % Magnitude of the final position vector
+
+ % Compute the cosine of delta nu directly from the position vectors
+ cos_delta_nu = dot(r0, rf) / (r0_mag * rf_mag); % Cosine of the angle between r0 and rf
+
+ % Initialize A to zero for the case where delta nu is zero
+ A = 0;
+ if cos_delta_nu < 1 % Only perform the calculation if the vectors are not parallel
+ if nargin < 5 || isempty(DM)
+ % Calculate delta nu using atan2 to ensure the correct quadrant
+ delta_nu = atan2(norm(cross(r0,rf)), dot(r0,rf));
+
+ % Compute DM based on delta nu
+ DM = 2 * (delta_nu < pi) - 1; % DM determines the direction of motion
+ end
+
+ % Calculate A using the formula given
+ A = DM * sqrt(r0_mag * rf_mag * (1 + cos_delta_nu)); % Scaled geometric parameter for the Lambert problem
+ end
+
+ % Initial guess for psi
+ psi = 0; % Initial guess for the universal anomaly related parameter
+ psi_low = -4 * pi^2; % Lower bound for psi
+ psi_up = 4 * pi^2; % Upper bound for psi
+ c2 = 0.5; % Initial guess for the C2 function, related to the universal variable
+ c3 = 1/6; % Initial guess for the C3 function, related to the universal variable
+
+ % Ensure we don't enter an infinite loop
+ while iter_count < max_iter % Loop until the solution converges or maximum iterations are reached
+ iter_count = iter_count + 1; % Increment iteration counter
+
+ [c2, c3] = stumpff(psi); % Calculate Stumpff functions C2 and C3 based on current psi
+ y = r0_mag + rf_mag - A * (1 - psi * c3) / sqrt(c2); % Calculate y, the universal variable related to the geometry of the trajectory
+
+ if A > 0 && y < 0 % Check if y is negative, which is physically infeasible
+ while y < 0 && iter_count < max_iter % Adjust psi upwards until y becomes positive or max_iter is reached
+ psi = psi + 0.1; % Increment psi by 0.1 to move toward a feasible region
+ [c2, c3] = stumpff(psi); % Recalculate Stumpff functions with new psi
+ y = r0_mag + rf_mag - A * (1 - psi * c3) / sqrt(c2); % Recalculate y with new psi
+ iter_count = iter_count + 1; % Increment iteration counter
+ end
+ end
+
+
+ chi = sqrt(y / c2); % Calculate the chi variable, a function of y and C2
+ Dt_old = (chi^3 * c3 + A * sqrt(y)) / sqrt(mu); % Compute the estimated time of flight based on current psi
+
+ if abs(Dt_old - Dt) < tol % Check if the calculated time of flight is within the tolerance of the desired time
+ exitflag = 1; % Set exitflag to 1 indicating successful convergence
+ break; % Exit the loop as convergence is achieved
+ elseif Dt_old < Dt % If the calculated time is less than desired, adjust psi lower bound
+ psi_low = psi;
+ else % If the calculated time is greater, adjust psi upper bound
+ psi_up = psi;
+ end
+ psi = (psi_up + psi_low) / 2; % Update psi to the average of upper and lower bounds for the next iteration
+ end
+
+ if iter_count >= max_iter
+ v0 = []; % Return empty arrays if solution was not found within the maximum iterations
+ vf = [];
+ exitflag = 0; % Set exitflag to 0 to indicate failure to converge
+ return;
+ end
+
+ % Calculate Lagrange coefficients f and g
+ f = 1 - y / r0_mag; % Calculate the f function, related to the geometry of the transfer orbit
+ g = A * sqrt(y / mu); % Calculate the g function, a scaling factor in the velocity equation
+ g_dot = 1 - y / rf_mag; % Calculate the derivative of g, another factor affecting velocity
+
+ % Calculate velocity vectors
+ v0 = (rf - f * r0) / g; % Calculate initial velocity vector required for the transfer
+ vf = (g_dot * rf - r0) / g; % Calculate final velocity vector at the end of the transfer
+
+
+ %% Stumpff Function
+ % Stumpff functions needed for the c2 and c3 calculations
+ function [c2, c3] = stumpff(psi)
+ if psi > tol % Check if psi is significantly greater than zero
+ % Use the series expansion for c2 when psi is positive
+ c2 = (1 - cos(sqrt(psi))) / psi;
+ % Calculate c3 using the sine series expansion for positive psi
+ c3 = (sqrt(psi) - sin(sqrt(psi))) / sqrt(psi^3);
+ elseif psi < -tol % Check if psi is significantly less than zero
+ % Use the hyperbolic cosine series for c2 when psi is negative
+ c2 = (1 - cosh(sqrt(-psi))) / psi;
+ % Calculate c3 using the hyperbolic sine series for negative psi
+ c3 = (sinh(sqrt(-psi)) - sqrt(-psi)) / sqrt((-psi)^3);
+ else % Handle the case where psi is approximately zero
+ % Default values based on the Taylor expansion around zero
+ c2 = 1/2;
+ c3 = 1/6;
+ end
+ end
+end
\ No newline at end of file