diff --git a/README.md b/README.md
index eb84a506..57a86ee1 100644
--- a/README.md
+++ b/README.md
@@ -4,18 +4,24 @@
 [![codecov](https://codecov.io/gh/leonlan/VRPLIB/branch/master/graph/badge.svg?token=X0X66LBNZ7)](https://codecov.io/gh/leonlan/VRPLIB)
 
 `vrplib` is a Python package for working with Vehicle Routing Problem (VRP) instances. The main features are:
-- reading VRPLIB and Solomon instances and solutions, and
+- reading VRPLIB and Solomon instances and solutions, 
+- writing VRPLIB-style instances and solutions, and
 - downloading instances and best known solutions from [CVRPLIB](http://vrp.atd-lab.inf.puc-rio.br/index.php/en/).
 
+## Outline
+- [Installation](#installation)
+- [Example usage](#example-usage)
+- [Documentation](#documentation)
+
 ## Installation
-`vrplib` works with Python 3.8+ and only depends on `numpy`.
+`vrplib` works with Python 3.8+ and only depends on `numpy`. It may be installed in the usual way as 
 
 ```shell
 pip install vrplib
 ```
 
 ## Example usage
-### Reading instances and solutions
+### Reading files
 ```python
 import vrplib
 
@@ -38,7 +44,72 @@ dict_keys(['routes', 'cost'])
 ```
 
 
-### Downloading instances from CVRPLIB 
+### Writing files
+The functions `write_instance` and `write_solution` provide a simple interface to writing instances and solutions in VRPLIB-style:
+- `write_instance` adds indices to data sections when necessary (`EDGE_WEIGHT_SECTION` and `DEPOT_SECTION` are excluded).
+- `write_solution` adds the `Route #{idx}` prefix to routes.
+
+Note that these functions do not validate instances: it is up to the user to write correct VRPLIB-style files.
+
+#### Instances
+``` python
+import vrplib
+
+instance_loc = "instance.vrp"
+instance_data = {
+    "NAME": "instance",
+    "TYPE": "CVRP",
+    "VEHICLES": 2,
+    "DIMENSION": 1,
+    "CAPACITY": 1,
+    "EDGE_WEIGHT_TYPE": "EUC_2D",
+    "NODE_COORD_SECTION": [[250, 250], [500, 500]],
+    "DEMAND_SECTION": [1, 1],
+    "DEPOT_SECTION": [1],
+}
+
+vrplib.write_instance(instance_loc, instance_data)
+```
+
+``` 
+NAME: instance
+TYPE: CVRP
+VEHICLES: 2
+DIMENSION: 1
+CAPACITY: 1
+EDGE_WEIGHT_TYPE: EUC_2D
+NODE_COORD_SECTION
+1	250	250
+2	500	500
+DEMAND_SECTION
+1	1
+2	1
+DEPOT_SECTION
+1
+EOF
+```
+
+#### Solutions
+``` python
+import vrplib
+
+solution_loc = "solution.sol"
+routes = [[1], [2, 3], [4, 5, 6]]
+solution_data = {"Cost": 42, "Vehicle types": [1, 2, 3]}
+
+vrplib.write_solution(solution_loc, routes, solution_data)
+```
+
+``` { .html } 
+Route #1: 1
+Route #2: 2 3
+Route #3: 4 5 6
+Cost: 42
+Vehicle types: [1, 2, 3]
+```
+
+
+### Downloading from CVRPLIB 
 ``` python
 import vrplib
 
diff --git a/tests/write/__init__.py b/tests/write/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/write/test_write_instance.py b/tests/write/test_write_instance.py
new file mode 100644
index 00000000..bac00606
--- /dev/null
+++ b/tests/write/test_write_instance.py
@@ -0,0 +1,157 @@
+import numpy as np
+from numpy.testing import assert_equal
+from pytest import mark
+
+from vrplib import write_instance
+
+
+@mark.parametrize(
+    "key, value, desired",
+    (
+        ["name", "Instance", "name: Instance"],  # string
+        ["DIMENSION", 100, "DIMENSION: 100"],  # int
+        ["VEHICLES", -10, "VEHICLES: -10"],  # negative
+        ["CAPACITY", 10.5, "CAPACITY: 10.5"],  # float
+        ["EMPTY", "", "EMPTY: "],  # empty
+    ),
+)
+def test_specifications(tmp_path, key, value, desired):
+    """
+    Tests that key-value pairs where values are floats or strings are
+    formatted as specifications.
+    """
+    name = "specifications"
+    instance = {key: value}
+    write_instance(tmp_path / name, instance)
+
+    desired = "\n".join([desired, "EOF", ""])
+    with open(tmp_path / name, "r") as fh:
+        assert_equal(fh.read(), desired)
+
+
+@mark.parametrize(
+    "key, value, desired",
+    (
+        # 1-dimensional list
+        ["X_SECTION", [0, 10], "\n".join(["X_SECTION", "1\t0", "2\t10"])],
+        # 1-dimensional list with mixed int and float values
+        ["X_SECTION", [0, 10.5], "\n".join(["X_SECTION", "1\t0", "2\t10.5"])],
+        # 1-dimensional list empty
+        ["X_SECTION", [], "\n".join(["X_SECTION"])],
+        # 2-dimensional numpy array
+        [
+            "Y_SECTION",
+            np.array([[0, 0], [1, 1]]),
+            "\n".join(["Y_SECTION", "1\t0\t0", "2\t1\t1"]),
+        ],
+        # 2-dimensional list empty
+        ["Y_SECTION", [[]], "\n".join(["Y_SECTION", "1\t"])],
+        # 2-dimensional array with different row lengths
+        # NOTE: This is currently an invalid VRPLIB format, see
+        # https://github.com/leonlan/VRPLIB/issues/108.
+        [
+            "DATA_SECTION",
+            [[1], [3, 4]],
+            "\n".join(["DATA_SECTION", "1\t1", "2\t3\t4"]),
+        ],
+    ),
+)
+def test_sections(tmp_path, key, value, desired):
+    """
+    Tests that key-value pairs where values are lists are formatted as
+    sections.
+    """
+    name = "sections"
+    instance = {key: value}
+    write_instance(tmp_path / name, instance)
+
+    with open(tmp_path / name, "r") as fh:
+        assert_equal(fh.read(), "\n".join([desired, "EOF", ""]))
+
+
+def test_no_indices_depot_and_edge_weight_section(tmp_path):
+    """
+    Tests that indices are not included when formatting depot and edge weight
+    section.
+    """
+    # Let's first test the depot section.
+    name = "depot"
+    instance = {"DEPOT_SECTION": [1, 2]}
+    write_instance(tmp_path / name, instance)
+
+    desired = "\n".join(["DEPOT_SECTION", "1", "2", "EOF", ""])
+    with open(tmp_path / name, "r") as fh:
+        assert_equal(fh.read(), desired)
+
+    # Now let's test the edge weight section.
+    name = "edge_weight"
+    instance = {
+        "EDGE_WEIGHT_SECTION": [
+            [1, 1, 2],
+            [1, 0, 3],
+            [1, 3, 0],
+        ]
+    }
+    write_instance(tmp_path / name, instance)
+
+    desired = "\n".join(
+        [
+            "EDGE_WEIGHT_SECTION",
+            "1\t1\t2",
+            "1\t0\t3",
+            "1\t3\t0",
+            "EOF",
+            "",
+        ]
+    )
+    with open(tmp_path / name, "r") as fh:
+        assert_equal(fh.read(), desired)
+
+
+def test_small_instance_example(tmp_path):
+    """
+    Tests if writing a small instance yields the correct result.
+    """
+    name = "C101"
+    instance = {
+        "NAME": name,
+        "TYPE": "VRPTW",
+        "DIMENSION": 4,
+        "CAPACITY": 200,
+        "NODE_COORD_SECTION": [
+            [40, 50],
+            [45, 68],
+            [45, 70],
+            [42, 66],
+        ],
+        "DEMAND_SECTION": [0, 10, 30, 10],
+        "DEPOT_SECTION": [1],
+    }
+
+    write_instance(tmp_path / name, instance)
+
+    desired = "\n".join(
+        [
+            "NAME: C101",
+            "TYPE: VRPTW",
+            "DIMENSION: 4",
+            "CAPACITY: 200",
+            "NODE_COORD_SECTION",
+            "1\t40\t50",
+            "2\t45\t68",
+            "3\t45\t70",
+            "4\t42\t66",
+            "DEMAND_SECTION",
+            "1\t0",
+            "2\t10",
+            "3\t30",
+            "4\t10",
+            "DEPOT_SECTION",
+            "1",
+            "EOF",
+            "",
+        ]
+    )
+
+    with open(tmp_path / name, "r") as fh:
+        assert_equal(fh.read(), desired)
diff --git a/tests/write/test_write_solution.py b/tests/write/test_write_solution.py
new file mode 100644
index 00000000..a5b65dd6
--- /dev/null
+++ b/tests/write/test_write_solution.py
@@ -0,0 +1,82 @@
+from numpy.testing import assert_equal, assert_raises
+from pytest import mark
+
+from vrplib import write_solution
+
+
+@mark.parametrize(
+    "routes, desired",
+    [
+        ([[1, 2]], "Route #1: 1 2"),
+        ([[1, 2], [42, 9]], "Route #1: 1 2\nRoute #2: 42 9"),
+    ],
+)
+def test_write_routes(tmp_path, routes, desired):
+    """
+    Tests the writing of a solution with routes.
+    """
+    name = "test.sol"
+    write_solution(tmp_path / name, routes)
+
+    with open(tmp_path / name, "r") as fh:
+        assert_equal(fh.read(), desired + "\n")
+
+
+def test_raise_empty_routes(tmp_path):
+    """
+    Tests that an error is raised if a route is empty.
+    """
+    name = "test.sol"
+
+    with assert_raises(ValueError):
+        write_solution(tmp_path / name, [[]])
+
+    with assert_raises(ValueError):
+        write_solution(tmp_path / name, [[1], []])
+
+
+@mark.parametrize(
+    "data, desired",
+    [
+        ({"Cost": 100}, "Cost: 100"),  # int
+        ({"Time": 123.45}, "Time: 123.45"),  # float
+        ({"Distance": -1}, "Distance: -1"),  # negative int
+        ({"name": "test.sol"}, "name: test.sol"),  # string
+        ({"Vehicle types": [1, 2, 3]}, "Vehicle types: [1, 2, 3]"),  # list
+        ({"Vehicle types": (1, 3)}, "Vehicle types: (1, 3)"),  # tuple
+    ],
+)
+def test_format_other_data(tmp_path, data, desired):
+    name = "test.sol"
+    routes = [[1]]
+    write_solution(tmp_path / name, routes, data)
+
+    with open(tmp_path / name, "r") as fh:
+        text = "Route #1: 1" + "\n" + desired + "\n"
+        assert_equal(fh.read(), text)
+
+
+def test_small_example(tmp_path):
+    """
+    Tests the writing of a small example.
+    """
+    name = "test.sol"
+    routes = [[1, 2], [3, 4], [5]]
+    data = {"Cost": 100, "Time": 123.45, "name": name}
+
+    write_solution(tmp_path / name, routes, data)
+
+    desired = "\n".join(
+        [
+            "Route #1: 1 2",
+            "Route #2: 3 4",
+            "Route #3: 5",
+            "Cost: 100",
+            "Time: 123.45",
+            "name: test.sol",
+            "",
+        ]
+    )
+
+    with open(tmp_path / name, "r") as fh:
+        assert_equal(fh.read(), desired)
diff --git a/vrplib/__init__.py b/vrplib/__init__.py
index 968e330b..380fe7ad 100644
--- a/vrplib/__init__.py
+++ b/vrplib/__init__.py
@@ -1,2 +1,3 @@
 from .download import download_instance, download_solution, list_names
 from .read import read_instance, read_solution
+from .write import write_instance, write_solution
diff --git a/vrplib/write/__init__.py b/vrplib/write/__init__.py
new file mode 100644
index 00000000..50a09a0b
--- /dev/null
+++ b/vrplib/write/__init__.py
@@ -0,0 +1,2 @@
+from .write_instance import write_instance
+from .write_solution import write_solution
diff --git a/vrplib/write/write_instance.py b/vrplib/write/write_instance.py
new file mode 100644
index 00000000..a022bdc6
--- /dev/null
+++ b/vrplib/write/write_instance.py
@@ -0,0 +1,94 @@
+import os
+from typing import Dict, List, Tuple, TypeVar, Union
+
+import numpy as np
+
+_ArrayLike = TypeVar("_ArrayLike", List, Tuple, np.ndarray)
+
+
+def write_instance(
+    path: Union[str, os.PathLike],
+    data: Dict[str, Union[str, int, float, _ArrayLike]],
+):
+    """
+    Writes a VRP instance to file following the VRPLIB format [1].
+
+    Parameters
+    ---------
+    path
+        The file path.
+    data
+        A dictionary of keyword-value pairs. For each key-value pair, the
+        following rules apply:
+        * If ``value`` is a string, integer or float, then it is considered a
+          problem specification and formatted as "{keyword}: {value}".
+        * If ``value`` is a one or two-dimensional array, then it is considered
+          a data section and formatted as
+          ```
+          {name}
+          1 {row_1}
+          2 {row_2}
+          ...
+          n {row_n}
+          ```
+          where ``name`` is the key and ``row_1``, ``row_2``, etc. are the
+          elements of the array. One-dimensional arrays are treated as column
+          vectors. If name is "EDGE_WEIGHT_SECTION" or "DEPOT_SECTION", then
+          the index is not included.
+
+    References
+    ----------
+    [1] Helsgaun, K. (2017). An Extension of the Lin-Kernighan-Helsgaun TSP
+        Solver for Constrained Traveling Salesman and Vehicle Routing
+        Problems.
+        http://webhotel4.ruc.dk/~keld/research/LKH-3/LKH-3_REPORT.pdf
+
+    """
+    with open(path, "w") as fh:
+        for key, value in data.items():
+            if isinstance(value, (str, int, float)):
+                fh.write(f"{key}: {value}" + "\n")
+            else:
+                fh.write(_format_section(key, value) + "\n")
+
+        fh.write("EOF\n")
+
+
+def _format_section(name: str, data: _ArrayLike) -> str:
+    """
+    Formats a data section.
+
+    Parameters
+    ----------
+    name
+        The name of the section.
+    data
+        The data to be formatted.
+
+    Returns
+    -------
+    str
+        A VRPLIB-formatted data section.
+    """
+    section = [name]
+    include_idx = name not in ["EDGE_WEIGHT_SECTION", "DEPOT_SECTION"]
+
+    if _is_one_dimensional(data):
+        # Treat 1D arrays as column vectors, so each element is a row.
+        for idx, elt in enumerate(data, 1):
+            prefix = f"{idx}\t" if include_idx else ""
+            section.append(prefix + str(elt))
+    else:
+        for idx, row in enumerate(data, 1):
+            prefix = f"{idx}\t" if include_idx else ""
+            rest = "\t".join([str(elt) for elt in row])
+            section.append(prefix + rest)
+
+    return "\n".join(section)
+
+
+def _is_one_dimensional(data: _ArrayLike) -> bool:
+    for elt in data:
+        if isinstance(elt, (list, tuple, np.ndarray)):
+            return False
+    return True
diff --git a/vrplib/write/write_solution.py b/vrplib/write/write_solution.py
new file mode 100644
index 00000000..75ee70b6
--- /dev/null
+++ b/vrplib/write/write_solution.py
@@ -0,0 +1,32 @@
+from typing import Any, Dict, List, Optional
+
+
+def write_solution(
+    path: str, routes: List[List[int]], data: Optional[Dict[str, Any]] = None
+):
+    """
+    Writes a VRP solution to file following the VRPLIB convention.
+
+    Parameters
+    ----------
+    path
+        The file path.
+    routes
+        A list of routes, each route denoting the order in which the customers
+        are visited.
+    **kwargs
+        Optional keyword arguments. Each key-value pair is written to the
+        solution file as "{key}: {value}".
+    """
+    for route in routes:
+        if len(route) == 0:
+            raise ValueError("Empty route in solution.")
+
+    with open(path, "w") as fi:
+        for idx, route in enumerate(routes, 1):
+            text = " ".join([f"Route #{idx}:"] + [str(val) for val in route])
+            fi.write(text + "\n")
+
+        if data is not None:
+            for key, value in data.items():
+                fi.write(f"{key}: {value}\n")