diff --git a/include/eal.h b/include/eal.h
index 423be9dcfa105290411a36976fe0b600047debfd..37bdd000b7bfd9bce4f740547c251e0dc6c0c9a3 100644
--- a/include/eal.h
+++ b/include/eal.h
@@ -38,7 +38,7 @@ namespace eal {
                                   std::vector<double> coefficients, double time);
   std::vector<double> getAngularVelocity(std::string from, std::string to,
                                          std::vector<double> coefficients, double time);
-  double evaluatePolynomial(std::vector<double> coeffs, double time);
+  double evaluatePolynomial(std::vector<double> coeffs, double time, int d);
   double interpolate(std::vector<double> points, std::vector<double> times, double time, interpolation interp, int d);
 }
 
diff --git a/src/eal.cpp b/src/eal.cpp
index 7d5f34cfa95963ccb871978b541fdd6887639657..490e5ec627fe2d3c55c6a73f9d17deab79fce4c1 100644
--- a/src/eal.cpp
+++ b/src/eal.cpp
@@ -85,9 +85,9 @@ namespace eal {
     }
 
     vector<double> coordinate = {0.0, 0.0, 0.0};
-    coordinate[0] = evaluatePolynomial(coeffs[0], time); // X
-    coordinate[1] = evaluatePolynomial(coeffs[1], time); // Y
-    coordinate[2] = evaluatePolynomial(coeffs[2], time); // Z
+    coordinate[0] = evaluatePolynomial(coeffs[0], time, 0); // X
+    coordinate[1] = evaluatePolynomial(coeffs[1], time, 0); // Y
+    coordinate[2] = evaluatePolynomial(coeffs[2], time, 0); // Z
 
     return coordinate;
   }
@@ -96,7 +96,16 @@ namespace eal {
   // Velocity Function
   // Takes the coefficients from the position equation
   vector<double> getVelocity(vector<vector<double>> coeffs, double time) {
+
+    if (coeffs.size() != 3) {
+      throw invalid_argument("Invalid input coeffs, expected three vectors.");
+    }
+
     vector<double> coordinate = {0.0, 0.0, 0.0};
+    coordinate[0] = evaluatePolynomial(coeffs[0], time, 1); // X
+    coordinate[1] = evaluatePolynomial(coeffs[1], time, 1); // Y
+    coordinate[2] = evaluatePolynomial(coeffs[2], time, 1); // Z
+
     return coordinate;
   }
 
@@ -186,13 +195,25 @@ namespace eal {
   // Polynomial evaluation helper function
   // The equation evaluated by this function is:
   //                x = cx_0 + cx_1 * t^(1) + ... + cx_n * t^n
-  double evaluatePolynomial(vector<double> coeffs, double time){
+  // The d parameter is for which derivative of the polynomial to compute.
+  // Supported options are
+  //   0: no derivative
+  //   1: first derivative
+  //   2: second derivative
+  double evaluatePolynomial(vector<double> coeffs, double time, int d){
     if (coeffs.empty()) {
       throw invalid_argument("Invalid input coeffs, must be non-empty.");
     }
 
-    const double *coeffsArray = coeffs.data();
-    return gsl_poly_eval(coeffsArray, coeffs.size(), time);
+    if (d < 0) {
+      throw invalid_argument("Invalid derivative degree, must be non-negative.");
+    }
+
+    vector<double> derivatives(d + 1);
+    gsl_poly_eval_derivs(coeffs.data(), coeffs.size(), time,
+                         derivatives.data(), derivatives.size());
+
+    return derivatives.back();
   }
 
  double interpolate(vector<double> points, vector<double> times, double time, interpolation interp, int d) {
diff --git a/tests/EalTest.cpp b/tests/EalTest.cpp
index 37be2dfe55e79928a61872f3394b38ef4ee53102..9f7c83f311c36634281684a678a5bb1bdc1a284b 100644
--- a/tests/EalTest.cpp
+++ b/tests/EalTest.cpp
@@ -136,6 +136,27 @@ TEST(SplineInterpTest, Extrapolate) {
                invalid_argument);
 }
 
+TEST(PolynomialTest, Evaluate) {
+  vector<double> coeffs = {1.0, 2.0, 3.0}; // 1 + 2x + 3x^2
+  EXPECT_EQ(2.0, eal::evaluatePolynomial(coeffs, -1, 0));
+}
+
+TEST(PolynomialTest, Derivatives) {
+  vector<double> coeffs = {1.0, 2.0, 3.0}; // 1 + 2x + 3x^2
+  EXPECT_EQ(-4.0, eal::evaluatePolynomial(coeffs, -1, 1));
+  EXPECT_EQ(6.0, eal::evaluatePolynomial(coeffs, -1, 2));
+}
+
+TEST(PolynomialTest, EmptyCoeffs) {
+  vector<double> coeffs = {};
+  EXPECT_THROW(eal::evaluatePolynomial(coeffs, -1, 1), invalid_argument);
+}
+
+TEST(PolynomialTest, BadDerivative) {
+  vector<double> coeffs = {1.0, 2.0, 3.0};
+  EXPECT_THROW(eal::evaluatePolynomial(coeffs, -1, -1), invalid_argument);
+}
+
 TEST(PoisitionCoeffTest, SecondOrderPolynomial) {
   double time = 2.0;
   vector<vector<double>> coeffs = {{1.0, 2.0, 3.0},
@@ -188,6 +209,30 @@ TEST(PoisitionCoeffTest, InvalidInput) {
 }
 
 
+TEST(VelocityCoeffTest, SecondOrderPolynomial) {
+  double time = 2.0;
+  vector<vector<double>> coeffs = {{1.0, 2.0, 3.0},
+                                   {1.0, 3.0, 2.0},
+                                   {3.0, 2.0, 1.0}};
+
+  vector<double> coordinate = eal::getVelocity(coeffs, time);
+
+  ASSERT_EQ(3, coordinate.size());
+  EXPECT_DOUBLE_EQ(14.0, coordinate[0]);
+  EXPECT_DOUBLE_EQ(11.0, coordinate[1]);
+  EXPECT_DOUBLE_EQ(6.0, coordinate[2]);
+}
+
+
+TEST(VelocityCoeffTest, InvalidInput) {
+  double valid_time = 0.0;
+  vector<vector<double>> invalid_coeffs_sizes = {{3.0, 2.0, 1.0},
+                                                 {1.0, 2.0, 3.0}};
+
+  EXPECT_THROW(eal::getVelocity(invalid_coeffs_sizes, valid_time), invalid_argument);
+}
+
+
 TEST(LinearInterpTest, ExmapleGetRotation) {
   // simple test, only checks if API hit correctly and output is normalized
   vector<double> times = {0,  1,  2, 3};