diff --git a/ext/or-tools/ext.cpp b/ext/or-tools/ext.cpp
index 2802836..eea9e4b 100644
--- a/ext/or-tools/ext.cpp
+++ b/ext/or-tools/ext.cpp
@@ -10,6 +10,7 @@ void init_assignment(Rice::Module& m);
 void init_bin_packing(Rice::Module& m);
 void init_constraint(Rice::Module& m);
 void init_linear(Rice::Module& m);
+void init_math_opt(Rice::Module& m);
 void init_network_flows(Rice::Module& m);
 void init_routing(Rice::Module& m);
 
@@ -24,6 +25,7 @@ void Init_ext()
   init_bin_packing(m);
   init_constraint(m);
   init_linear(m);
+  init_math_opt(m);
   init_network_flows(m);
   init_routing(m);
 
diff --git a/ext/or-tools/math_opt.cpp b/ext/or-tools/math_opt.cpp
new file mode 100644
index 0000000..ace71e1
--- /dev/null
+++ b/ext/or-tools/math_opt.cpp
@@ -0,0 +1,25 @@
+#include "absl/log/check.h"
+#include "absl/status/statusor.h"
+#include "ortools/base/init_google.h"
+#include "ortools/math_opt/cpp/math_opt.h"
+
+#include "ext.h"
+
+using operations_research::math_opt::BoundedLinearExpression;
+using operations_research::math_opt::Model;
+using operations_research::math_opt::Variable;
+
+void init_math_opt(Rice::Module& m) {
+  auto mathopt = Rice::define_module_under(m, "MathOpt");
+
+  Rice::define_class_under<Variable>(mathopt, "Variable");
+
+  Rice::define_class_under<Model>(mathopt, "Model")
+    .define_constructor(Rice::Constructor<Model, std::string>())
+    .define_method("add_variable", &Model::AddContinuousVariable)
+    .define_method(
+      "add_linear_constraint",
+      [](Model& self, const BoundedLinearExpression& bounded_expr, const std::string& name) {
+        self.AddLinearConstraint(bounded_expr, name);
+      });
+}
diff --git a/test/math_opt_test.rb b/test/math_opt_test.rb
new file mode 100644
index 0000000..8db4dcf
--- /dev/null
+++ b/test/math_opt_test.rb
@@ -0,0 +1,27 @@
+require_relative "test_helper"
+
+class MathOptTest < Minitest::Test
+  # https://developers.google.com/optimization/math_opt/basic_example
+  def test_basic
+    model = ORTools::MathOpt::Model.new("getting_started_lp")
+    x = model.add_variable(-1.0, 1.5, "x")
+    y = model.add_variable(0.0, 1.0, "y")
+
+    skip "todo"
+
+    model.add_linear_constraint(x + y <= 1.5)
+    model.maximize(x + 2 * y)
+
+    params = ORTools::MathOpt::SolveParameters.new(enable_output: true)
+
+    result = mathopt.solve(model, :glop, params)
+    if result.termination.reason != :optimal
+      raise RuntimeError, "model failed to solve: #{result.termination}"
+    end
+
+    puts "MathOpt solve succeeded"
+    puts "Objective value:", result.objective_value
+    puts "x:", result.variable_values[x]
+    puts "y:", result.variable_values[y]
+  end
+end
diff --git a/test/support/mathopt.py b/test/support/mathopt.py
new file mode 100644
index 0000000..a1fc57a
--- /dev/null
+++ b/test/support/mathopt.py
@@ -0,0 +1,24 @@
+from ortools.math_opt.python import mathopt
+
+# Build the model.
+model = mathopt.Model(name="getting_started_lp")
+x = model.add_variable(lb=-1.0, ub=1.5, name="x")
+y = model.add_variable(lb=0.0, ub=1.0, name="y")
+model.add_linear_constraint(x + y <= 1.5)
+model.maximize(x + 2 * y)
+
+# Set parameters, e.g. turn on logging.
+params = mathopt.SolveParameters(enable_output=True)
+
+# Solve and ensure an optimal solution was found with no errors.
+# (mathopt.solve may raise a RuntimeError on invalid input or internal solver
+# errors.)
+result = mathopt.solve(model, mathopt.SolverType.GLOP, params=params)
+if result.termination.reason != mathopt.TerminationReason.OPTIMAL:
+    raise RuntimeError(f"model failed to solve: {result.termination}")
+
+# Print some information from the result.
+print("MathOpt solve succeeded")
+print("Objective value:", result.objective_value())
+print("x:", result.variable_values()[x])
+print("y:", result.variable_values()[y])