diff --git a/numerical_methods/ode_forward_euler.cpp b/numerical_methods/ode_forward_euler.cpp
new file mode 100644
index 000000000..e289d96e7
--- /dev/null
+++ b/numerical_methods/ode_forward_euler.cpp
@@ -0,0 +1,207 @@
+/**
+ * \file
+ * \authors [Krishna Vedala](https://github.com/kvedala)
+ * \brief Solve a multivariable first order [ordinary differential equation
+ * (ODEs)](https://en.wikipedia.org/wiki/Ordinary_differential_equation) using
+ * [forward Euler
+ * method](https://en.wikipedia.org/wiki/Numerical_methods_for_ordinary_differential_equations#Euler_method)
+ *
+ * \details
+ * The ODE being solved is:
+ * \f{eqnarray*}{
+ * \dot{u} &=& v\\
+ * \dot{v} &=& -\omega^2 u\\
+ * \omega &=& 1\\
+ * [x_0, u_0, v_0] &=& [0,1,0]\qquad\ldots\text{(initial values)}
+ * \f}
+ * The exact solution for the above problem is:
+ * \f{eqnarray*}{
+ * u(x) &=& \cos(x)\\
+ * v(x) &=& -\sin(x)\\
+ * \f}
+ * The computation results are stored to a text file `forward_euler.csv` and the
+ * exact soltuion results in `exact.csv` for comparison.
+ * 
+ *
+ * To implement [Van der Pol
+ * oscillator](https://en.wikipedia.org/wiki/Van_der_Pol_oscillator), change the
+ * ::problem function to:
+ * ```cpp
+ * const double mu = 2.0;
+ * dy[0] = y[1];
+ * dy[1] = mu * (1.f - y[0] * y[0]) * y[1] - y[0];
+ * ```
+ * \see ode_midpoint_euler.cpp, ode_semi_implicit_euler.cpp
+ */
+
+#include
+#include
+#include
+#include
+#include
+
+/**
+ * @brief Problem statement for a system with first-order differential
+ * equations. Updates the system differential variables.
+ * \note This function can be updated to and ode of any order.
+ *
+ * @param[in] x independent variable(s)
+ * @param[in,out] y dependent variable(s)
+ * @param[in,out] dy first-derivative of dependent variable(s)
+ */
+void problem(const double &x, std::valarray *y,
+ std::valarray *dy) {
+ const double omega = 1.F; // some const for the problem
+ dy[0][0] = y[0][1]; // x dot
+ dy[0][1] = -omega * omega * y[0][0]; // y dot
+}
+
+/**
+ * @brief Exact solution of the problem. Used for solution comparison.
+ *
+ * @param[in] x independent variable
+ * @param[in,out] y dependent variable
+ */
+void exact_solution(const double &x, std::valarray *y) {
+ y[0][0] = std::cos(x);
+ y[0][1] = -std::sin(x);
+}
+
+/** \addtogroup ode Ordinary Differential Equations
+ * @{
+ */
+/**
+ * @brief Compute next step approximation using the forward-Euler
+ * method. @f[y_{n+1}=y_n + dx\cdot f\left(x_n,y_n\right)@f]
+ * @param[in] dx step size
+ * @param[in] x take \f$x_n\f$ and compute \f$x_{n+1}\f$
+ * @param[in,out] y take \f$y_n\f$ and compute \f$y_{n+1}\f$
+ * @param[in,out] dy compute \f$f\left(x_n,y_n\right)\f$
+ */
+void forward_euler_step(const double dx, const double &x,
+ std::valarray *y, std::valarray *dy) {
+ problem(x, y, dy);
+ y[0] += dy[0] * dx;
+}
+
+/**
+ * @brief Compute approximation using the forward-Euler
+ * method in the given limits.
+ * @param[in] dx step size
+ * @param[in] x0 initial value of independent variable
+ * @param[in] x_max final value of independent variable
+ * @param[in,out] y take \f$y_n\f$ and compute \f$y_{n+1}\f$
+ * @param[in] save_to_file flag to save results to a CSV file (1) or not (0)
+ * @returns time taken for computation in seconds
+ */
+double forward_euler(double dx, double x0, double x_max,
+ std::valarray *y, bool save_to_file = false) {
+ std::valarray dy = y[0];
+
+ std::ofstream fp;
+ if (save_to_file) {
+ fp.open("forward_euler.csv", std::ofstream::out);
+ if (!fp.is_open()) {
+ std::perror("Error! ");
+ }
+ }
+
+ std::size_t L = y->size();
+
+ /* start integration */
+ std::clock_t t1 = std::clock();
+ double x = x0;
+ do // iterate for each step of independent variable
+ {
+ if (save_to_file && fp.is_open()) {
+ // write to file
+ fp << x << ",";
+ for (int i = 0; i < L - 1; i++) {
+ fp << y[0][i] << ",";
+ }
+ fp << y[0][L - 1] << "\n";
+ }
+
+ forward_euler_step(dx, x, y, &dy); // perform integration
+ x += dx; // update step
+ } while (x <= x_max); // till upper limit of independent variable
+ /* end of integration */
+ std::clock_t t2 = std::clock();
+
+ if (fp.is_open())
+ fp.close();
+
+ return static_cast(t2 - t1) / CLOCKS_PER_SEC;
+}
+
+/** @} */
+
+/**
+ * Function to compute and save exact solution for comparison
+ *
+ * \param [in] X0 initial value of independent variable
+ * \param [in] X_MAX final value of independent variable
+ * \param [in] step_size independent variable step size
+ * \param [in] Y0 initial values of dependent variables
+ */
+void save_exact_solution(const double &X0, const double &X_MAX,
+ const double &step_size,
+ const std::valarray &Y0) {
+ double x = X0;
+ std::valarray y = Y0;
+
+ std::ofstream fp("exact.csv", std::ostream::out);
+ if (!fp.is_open()) {
+ std::perror("Error! ");
+ return;
+ }
+ std::cout << "Finding exact solution\n";
+
+ std::clock_t t1 = std::clock();
+ do {
+ fp << x << ",";
+ for (int i = 0; i < y.size() - 1; i++) {
+ fp << y[i] << ",";
+ }
+ fp << y[y.size() - 1] << "\n";
+
+ exact_solution(x, &y);
+
+ x += step_size;
+ } while (x <= X_MAX);
+
+ std::clock_t t2 = std::clock();
+ double total_time = static_cast(t2 - t1) / CLOCKS_PER_SEC;
+ std::cout << "\tTime = " << total_time << " ms\n";
+
+ fp.close();
+}
+
+/**
+ * Main Function
+ */
+int main(int argc, char *argv[]) {
+ double X0 = 0.f; /* initial value of x0 */
+ double X_MAX = 10.F; /* upper limit of integration */
+ std::valarray Y0 = {1.f, 0.f}; /* initial value Y = y(x = x_0) */
+ double step_size;
+
+ if (argc == 1) {
+ std::cout << "\nEnter the step size: ";
+ std::cin >> step_size;
+ } else {
+ // use commandline argument as independent variable step size
+ step_size = std::atof(argv[1]);
+ }
+
+ // get approximate solution
+ double total_time = forward_euler(step_size, X0, X_MAX, &Y0, true);
+ std::cout << "\tTime = " << total_time << " ms\n";
+
+ /* compute exact solution for comparion */
+ save_exact_solution(X0, X_MAX, step_size, Y0);
+
+ return 0;
+}