diff --git a/src/simulation/core_sim.zig b/src/simulation/core_sim.zig
new file mode 100644
index 0000000..2ca2068
--- /dev/null
+++ b/src/simulation/core_sim.zig
@@ -0,0 +1,24 @@
+const std = @import("std");
+const ODE = @import("../solvers/ode.zig.").ODE;
+const IntegrationMethod = @import("../solvers/ode.zig.").IntegrationMethod;
+const SolverError = @import("../solvers/ode.zig.").SolverError;
+
+pub fn simulate(StateDimension: usize, InputDimension: usize, initial_dt: f64, end_time: f64, derivative: fn (@Vector(StateDimension, f64)) @Vector(StateDimension, f64), method: IntegrationMethod) !void {
+    // Create initial state (assuming 0s for now)
+    var state = ODE(StateDimension, InputDimension).init(
+        @Vector([]f64, StateDimension),
+        null,
+        initial_dt,
+        0.0,
+        derivative,
+    );
+
+    var current_time = 0.0;
+    while (current_time < end_time) {
+        state = try state.step(null, method, null, null, null);
+
+        std.debug.print("Time: {}, State: {}\n", .{ state.t, state.x });
+
+        current_time += state.dt;
+    }
+}
diff --git a/src/simulation/test_sim.zig b/src/simulation/test_sim.zig
new file mode 100644
index 0000000..ab1e530
--- /dev/null
+++ b/src/simulation/test_sim.zig
@@ -0,0 +1,28 @@
+const std = @import("std");
+const core_sim = @import("core_sim.zig");
+const ODE = @import("../solvers/ode.zig.").ODE;
+const IntegrationMethod = @import("../solvers/ode.zig.").IntegrationMethod;
+
+const StateDimension = 2;
+const InputDimension = 1;
+
+fn simple_derivative(x: @Vector(StateDimension, f64)) @Vector(StateDimension, f64) {
+    // Exampe derivative: dx/dt = -x
+    var result: @Vector(StateDimension, f64) = undefined;
+    result[0] = -x[0];
+    result[1] = -x[1];
+    return result;
+}
+
+pub fn test_simulation() void {
+    const initial_dt = 0.1;
+    const end_time = 2.0;
+
+    // using Euler method
+    std.debug.print("Running simulation with Euler method...\n", .{});
+    _ = core_sim.simulate(StateDimension, InputDimension, initial_dt, end_time, simple_derivative, IntegrationMethod.euler);
+
+    // using Runge-Kutta method
+    std.debug.print("Running simulation with Runge-Kutta method...\n", .{});
+    _ = core_sim.simulate(StateDimension, InputDimension, initial_dt, end_time, simple_derivative, IntegrationMethod.rk4);
+}