diff --git a/sequential.cpp b/sequential.cpp
new file mode 100644
index 0000000..611cb56
--- /dev/null
+++ b/sequential.cpp
@@ -0,0 +1,106 @@
+#include <iostream>
+#include <vector>
+#include <iomanip>
+#include <math.h>
+
+using namespace std;
+
+int main() 
+{
+    // Length
+    int l = 1925;
+    // Spatial step
+    int dx = 40;
+    // Temporal step
+    int dt = 2500;
+    // Number of spatial steps
+    int kmax = l / dx;
+
+    // Initializing parameters
+    vector<double> list(dt);
+
+    // Arrays to store data
+    double Q[kmax][dt];
+    double P[kmax][dt];
+
+    int i;
+
+    // Parameters
+    double alpha = 0.093;
+    double beta = 0.004;
+    double gamma = 728.6;
+    double h = 0.05;
+    double sum_of_h = 0;
+
+    // Loop over time steps
+    for (i = 0; i < dt; i++) 
+    {
+        // Loop over spatial steps
+        for (int k = 0; k < kmax; k++) 
+        {
+            // Boundary condition for the last spatial step
+            if (k == kmax - 1) 
+            {
+                Q[k][i] = 0;
+                P[k][i] = -202.0; // P = -202 when Q = 10
+            } else 
+            {
+                Q[k][i] = 0;
+                P[k][i] = 0;
+            }
+        }
+    }
+
+    // Main calculation loop
+    for (i = 2; i < dt - 2; i++) 
+    {
+        // Add current time step to the list
+        sum_of_h = sum_of_h + h;
+        list.push_back(sum_of_h);
+
+		std::cout << "Q at step " << i;
+		
+        // Loop over spatial steps
+        for (int k = 1; k < kmax; k++) 
+        {
+            // Boundary condition for the last spatial step
+            if (k == kmax - 1) 
+            {
+                Q[k][i + 1] = Q[k][i] + h * (alpha * (P[k - 1][i] - P[k][i]) - beta * Q[k][i] * abs(Q[k][i]));
+            } else 
+            {
+                // Calculate Q and P for non-boundary spatial steps
+                Q[k][i + 1] = Q[k][i] + h * (alpha * (P[k - 1][i] - P[k][i]) - beta * Q[k][i] * abs(Q[k][i]));
+                P[k][i + 1] = P[k][i - 2] + h * (gamma * (Q[k][i] - Q[k + 1][i]));
+            }
+        }
+    }
+    
+    // Print out the first 10 values of the timestemp
+    std::cout << "FIRST 10 VALUES\n";
+    for(int i = 0; i < kmax; i++)
+    {
+        std::cout.width(5);
+        std::cout.flags(std::ios::left);
+        std::cout << i << ": ";
+   
+        for(int k = 0; k < 10; k++)
+            std::cout << std::setw(10) << Q[i][k];
+        std::cout << std::endl;
+    }
+
+    // Print out the last 10 values of the timestemp
+    std::cout << "LAST 10 VALUES\n";
+    for(int i = 0; i < kmax; i++)
+    {
+        std::cout.width(5);
+        std::cout.flags(std::ios::left);
+        std::cout << i << ": ";
+ 
+        for(int k = 0; k < 10; k++)
+            std::cout << std::setw(10) << Q[i][2489 + k];
+        std::cout << std::endl;
+    }
+
+    return 0;
+}