diff --git a/R/model.R b/R/model.R
index f9b8241c..0a1cb1c1 100644
--- a/R/model.R
+++ b/R/model.R
@@ -54,6 +54,9 @@
 #' file <- file.path(cmdstan_path(), "examples/bernoulli/bernoulli.stan")
 #' mod <- cmdstan_model(file)
 #' mod$print()
+#' # Print with line numbers. This can be set globally using the
+#' # `cmdstanr_print_line_numbers` option.
+#' mod$print(line_numbers = TRUE)
 #'
 #' # Data as a named list (like RStan)
 #' stan_data <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))
@@ -298,11 +301,21 @@ CmdStanModel <- R6::R6Class(
       }
       private$stan_code_
     },
-    print = function() {
+    print = function(line_numbers = getOption("cmdstanr_print_line_numbers", FALSE)) {
       if (length(private$stan_code_) == 0) {
         stop("'$print()' cannot be used because the 'CmdStanModel' was not created with a Stan file.", call. = FALSE)
       }
-      cat(self$code(), sep = "\n")
+      lines <- self$code()
+      if (line_numbers) {
+        line_num_indent <- nchar(as.character(length(lines)))
+        line_nums <- vapply(seq_along(lines), function(y) {
+          paste0(
+            rep(" ", line_num_indent - nchar(as.character(y))), y, collapse = ""
+          )
+        }, character(1))
+        lines <- paste(paste(line_nums, lines, sep = ": "), collapse = "\n")
+      }
+      cat(lines, sep = "\n")
       invisible(self)
     },
     stan_file = function() {
diff --git a/man/CmdStanModel.Rd b/man/CmdStanModel.Rd
index 55200357..8bb3fdaf 100644
--- a/man/CmdStanModel.Rd
+++ b/man/CmdStanModel.Rd
@@ -73,6 +73,9 @@ set_cmdstan_path(path = NULL)
 file <- file.path(cmdstan_path(), "examples/bernoulli/bernoulli.stan")
 mod <- cmdstan_model(file)
 mod$print()
+# Print with line numbers. This can be set globally using the
+# `cmdstanr_print_line_numbers` option.
+mod$print(line_numbers = TRUE)
 
 # Data as a named list (like RStan)
 stan_data <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))
diff --git a/man/cmdstan_model.Rd b/man/cmdstan_model.Rd
index b0432abc..bca22a30 100644
--- a/man/cmdstan_model.Rd
+++ b/man/cmdstan_model.Rd
@@ -59,6 +59,9 @@ set_cmdstan_path(path = NULL)
 file <- file.path(cmdstan_path(), "examples/bernoulli/bernoulli.stan")
 mod <- cmdstan_model(file)
 mod$print()
+# Print with line numbers. This can be set globally using the
+# `cmdstanr_print_line_numbers` option.
+mod$print(line_numbers = TRUE)
 
 # Data as a named list (like RStan)
 stan_data <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))
diff --git a/man/cmdstanr-package.Rd b/man/cmdstanr-package.Rd
index 739213ff..eb8ae80c 100644
--- a/man/cmdstanr-package.Rd
+++ b/man/cmdstanr-package.Rd
@@ -89,6 +89,9 @@ set_cmdstan_path(path = NULL)
 file <- file.path(cmdstan_path(), "examples/bernoulli/bernoulli.stan")
 mod <- cmdstan_model(file)
 mod$print()
+# Print with line numbers. This can be set globally using the
+# `cmdstanr_print_line_numbers` option.
+mod$print(line_numbers = TRUE)
 
 # Data as a named list (like RStan)
 stan_data <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))
diff --git a/man/model-method-optimize.Rd b/man/model-method-optimize.Rd
index 55d7758f..ac8f9f65 100644
--- a/man/model-method-optimize.Rd
+++ b/man/model-method-optimize.Rd
@@ -245,6 +245,9 @@ set_cmdstan_path(path = NULL)
 file <- file.path(cmdstan_path(), "examples/bernoulli/bernoulli.stan")
 mod <- cmdstan_model(file)
 mod$print()
+# Print with line numbers. This can be set globally using the
+# `cmdstanr_print_line_numbers` option.
+mod$print(line_numbers = TRUE)
 
 # Data as a named list (like RStan)
 stan_data <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))
diff --git a/man/model-method-pathfinder.Rd b/man/model-method-pathfinder.Rd
index c3a2044b..90578282 100644
--- a/man/model-method-pathfinder.Rd
+++ b/man/model-method-pathfinder.Rd
@@ -270,6 +270,9 @@ set_cmdstan_path(path = NULL)
 file <- file.path(cmdstan_path(), "examples/bernoulli/bernoulli.stan")
 mod <- cmdstan_model(file)
 mod$print()
+# Print with line numbers. This can be set globally using the
+# `cmdstanr_print_line_numbers` option.
+mod$print(line_numbers = TRUE)
 
 # Data as a named list (like RStan)
 stan_data <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))
diff --git a/man/model-method-sample.Rd b/man/model-method-sample.Rd
index bc5b2a8c..9b0b1ac4 100644
--- a/man/model-method-sample.Rd
+++ b/man/model-method-sample.Rd
@@ -348,6 +348,9 @@ set_cmdstan_path(path = NULL)
 file <- file.path(cmdstan_path(), "examples/bernoulli/bernoulli.stan")
 mod <- cmdstan_model(file)
 mod$print()
+# Print with line numbers. This can be set globally using the
+# `cmdstanr_print_line_numbers` option.
+mod$print(line_numbers = TRUE)
 
 # Data as a named list (like RStan)
 stan_data <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))
diff --git a/man/model-method-variational.Rd b/man/model-method-variational.Rd
index 3892d886..417a647a 100644
--- a/man/model-method-variational.Rd
+++ b/man/model-method-variational.Rd
@@ -245,6 +245,9 @@ set_cmdstan_path(path = NULL)
 file <- file.path(cmdstan_path(), "examples/bernoulli/bernoulli.stan")
 mod <- cmdstan_model(file)
 mod$print()
+# Print with line numbers. This can be set globally using the
+# `cmdstanr_print_line_numbers` option.
+mod$print(line_numbers = TRUE)
 
 # Data as a named list (like RStan)
 stan_data <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))