diff --git a/tests/hello481/breakpoints.seer b/tests/hello481/breakpoints.seer
index 50ee722..618f56b 100644
--- a/tests/hello481/breakpoints.seer
+++ b/tests/hello481/breakpoints.seer
@@ -1 +1 @@
-break -source /nas/erniep/Development/seer/tests/hello481/hello481.cpp -line 17
+break -source /nas/erniep/Development/seer/tests/hello481/hello481.cpp -line 81
diff --git a/tests/hello481/hello481.cpp b/tests/hello481/hello481.cpp
index 10bcfc4..bb6ebd8 100644
--- a/tests/hello481/hello481.cpp
+++ b/tests/hello481/hello481.cpp
@@ -1,19 +1,92 @@
-#include <string>
+// damped_oscillator
+//
+// Simulates a damped harmonic oscillator:
+//     x'' + 2*zeta*omega0*x' + omega0^2*x = 0
+//
+// Stores a 4000-point dataset directly into two double arrays:
+//     X[i] = position at step i
+//     Y[i] = velocity at step i
+//
+
 #include <iostream>
+#include <cmath>
+#include <iomanip>
 
-int main () {
+const int N = 4000;   // number of data points
 
-    double pos[40];
-    double vit[40];
+double X[N];  // position array
+double Y[N];  // velocity array
 
-    double posv = 1.0;
+struct State {
+    double x;  // position
+    double v;  // velocity
+};
 
-    for (int i=0; i<sizeof(pos); i++) {
+class DampedOscillator {
+    public:
+        DampedOscillator(double omega0, double zeta) : omega0_(omega0), zeta_(zeta) {}
 
-        pos[i] = posv;
-        posv -= (1.0 - 0.0) / double(sizeof(pos));
+        // dx/dt, dv/dt
+        State derivatives(const State& s) const {
+            State d;
+            d.x = s.v;
+            d.v = -2.0 * zeta_ * omega0_ * s.v - omega0_ * omega0_ * s.x;
+            return d;
+        }
+
+        // Single RK4 step
+        State step(const State& s, double dt) const {
+            State k1 = derivatives(s);
+
+            State s2{ s.x + 0.5 * dt * k1.x, s.v + 0.5 * dt * k1.v };
+            State k2 = derivatives(s2);
+
+            State s3{ s.x + 0.5 * dt * k2.x, s.v + 0.5 * dt * k2.v };
+            State k3 = derivatives(s3);
+
+            State s4{ s.x + dt * k3.x, s.v + dt * k3.v };
+            State k4 = derivatives(s4);
+
+            State result;
+            result.x = s.x + (dt / 6.0) * (k1.x + 2.0 * k2.x + 2.0 * k3.x + k4.x);
+            result.v = s.v + (dt / 6.0) * (k1.v + 2.0 * k2.v + 2.0 * k3.v + k4.v);
+            return result;
+        }
+
+    private:
+        double omega0_; // natural angular frequency (rad/s)
+        double zeta_;   // damping ratio (0 = undamped, 1 = critically damped)
+};
+
+int main() {
+
+    // ---- Simulation parameters (tweak as you like) ----
+    const double dt      = 0.01;   // time step (s)
+    const double omega0  = 2.0;    // natural angular frequency (rad/s)
+    const double zeta    = 0.05;   // damping ratio (underdamped)
+    const double x0      = 1.0;    // initial position
+    const double v0      = 0.0;    // initial velocity
+
+    DampedOscillator osc(omega0, zeta);
+    State state{ x0, v0 };
+
+    // Fill the arrays instead of writing to a file
+    for (int i = 0; i < N; ++i) {
+        X[i] = state.x;
+        Y[i] = state.v;
+        state = osc.step(state, dt);
     }
 
+    // Optional: print a few samples to confirm it worked
+    std::cout << std::fixed << std::setprecision(6);
+    std::cout << "First 5 points (X = position, Y = velocity):\n";
+
+    for (int i = 0; i < 5; ++i) {
+        std::cout << "  X[" << i << "]=" << X[i] << "  Y[" << i << "]=" << Y[i] << "\n";
+    }
+
+    std::cout << "Last point: X[" << N-1 << "]=" << X[N-1] << "  Y[" << N-1 << "]=" << Y[N-1] << "\n";
+
     return 0;
 }