diff --git a/.github/workflows/test_pytest.yaml b/.github/workflows/test_pytest.yaml
index d8bd00f70..ebd84b465 100644
--- a/.github/workflows/test_pytest.yaml
+++ b/.github/workflows/test_pytest.yaml
@@ -13,6 +13,7 @@ jobs:
           - windows-latest
         python-version:
           - 3.7
+          - 3.11
     steps:
     - uses: actions/checkout@v2
     - name: Set up Python ${{ matrix.python-version }}
diff --git a/docs/conf.py b/docs/conf.py
index 44d9a48d9..925ea26d1 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -24,7 +24,7 @@
 author = "Giovani Hidalgo Ceotto"
 
 # The full version, including alpha/beta/rc tags
-release = "0.13.0"
+release = "0.13.1"
 
 
 # -- General configuration ---------------------------------------------------
diff --git a/docs/user/installation.rst b/docs/user/installation.rst
index 3a3640930..fb9a64d2d 100644
--- a/docs/user/installation.rst
+++ b/docs/user/installation.rst
@@ -19,7 +19,7 @@ If you want to choose a specific version to guarantee compatibility, you may ins
 
 .. code-block:: shell
 
-    pip install rocketpy==0.13.0
+    pip install rocketpy==0.13.1
 
 
 Optional Installation Method: ``conda``
diff --git a/docs/user/requirements.rst b/docs/user/requirements.rst
index fe7b9f1c3..a6a5c07ef 100644
--- a/docs/user/requirements.rst
+++ b/docs/user/requirements.rst
@@ -8,7 +8,6 @@ Python Version
 --------------
 
 RocketPy supports Python 3.7 and above.
-Support for Python 3.11 is still limited by some dependencies.
 Sorry, there are currently no plans to support earlier versions.
 If you really need to run RocketPy on Python 3.6 or earlier, feel free to submit an issue and we will see what we can do!
 
@@ -21,7 +20,7 @@ The following packages are needed in order to run RocketPy:
 - Numpy >= 1.0
 - Scipy >= 1.0
 - Matplotlib >= 3.0
-- netCDF4 >= 1.4
+- netCDF4 >= 1.4, < 1.6 for Python 3.7+, netCDF4 >= 1.6.2 for Python 3.11
 - windrose >= 1.6.8
 - requests
 - pytz
@@ -29,10 +28,11 @@ The following packages are needed in order to run RocketPy:
 - ipywidgets >= 7.6.3
 - jsonpickle
 
- 
 All of these packages, are automatically installed when RocketPy is installed using either ``pip`` or ``conda``.
 However, in case the user wants to install these packages manually, they can do so by following the instructions bellow.
 
+Note: Google Colab and netCDF4 don't play well together for netCDF4 >= 1.6 (see `this issue <https://github.com/Unidata/netcdf4-python/issues/1179>`_). If you are using Google Colab, you should install netCDF4 < 1.6 or use `condacolab <https://github.com/conda-incubator/condacolab>`_) to install netCDF4 >= 1.6.
+
 Installing Required Packages Using ``pip``
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -43,7 +43,7 @@ The packages needed can be installed via ``pip`` by running the following lines
     pip install "numpy>=1.0" 
     pip install "scipy>=1.0"
     pip install "matplotlib>=3.0"
-    pip install "netCDF4>=1.4"
+    pip install "netCDF4>=1.6.2"
     pip install "windrose >= 1.6.8"
     pip install "ipywidgets>=7.6.3"
     pip install requests
@@ -61,7 +61,7 @@ To update Scipy and install netCDF4 using Conda, the following code is used:
 .. code-block:: shell
 
     conda install "scipy>=1.0"
-    conda install -c anaconda "netcdf4>=1.4"
+    conda install -c anaconda "netcdf4>=1.6.2"
 
 
 Optional Packages
diff --git a/requirements.txt b/requirements.txt
index ccfdb4f15..3d9fc3e49 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,7 +1,8 @@
 numpy>=1.0
 scipy>=1.0
 matplotlib>=3.0
-netCDF4>=1.4
+netCDF4>=1.4,<1.6; python_version == "3.7"
+netCDF4>=1.6.2; python_version != "3.7"
 windrose>=1.6.8
 ipywidgets>=7.6.3
 requests
diff --git a/rocketpy/AeroSurfaces.py b/rocketpy/AeroSurfaces.py
index fcd3064a9..e99562551 100644
--- a/rocketpy/AeroSurfaces.py
+++ b/rocketpy/AeroSurfaces.py
@@ -1414,7 +1414,6 @@ def geometricInfo(self):
         return None
 
     def aerodynamicInfo(self):
-
         print(f"\nTail name: {self.name}")
         print(f"Tail Center of Pressure: {self.cp}")
         print(f"Tail Lift Coefficient Slope: {self.clalpha}")
diff --git a/rocketpy/EnvironmentAnalysis.py b/rocketpy/EnvironmentAnalysis.py
index e318e043c..85b4d8327 100644
--- a/rocketpy/EnvironmentAnalysis.py
+++ b/rocketpy/EnvironmentAnalysis.py
@@ -637,7 +637,6 @@ def parsePressureLevelData(self):
 
         # Loop through time and save all values
         for timeIndex, timeNum in enumerate(timeNumArray):
-
             dateString, hourString, dateTime = self.__timeNumToDateString(
                 timeNum, timeNumArray.units, calendar="gregorian"
             )
@@ -820,7 +819,6 @@ def parseSurfaceData(self):
 
         # Loop through time and save all values
         for timeIndex, timeNum in enumerate(timeNumArray):
-
             dateString, hourString, dateTime = self.__timeNumToDateString(
                 timeNum, timeNumArray.units, calendar="gregorian"
             )
@@ -2130,7 +2128,7 @@ def plot_wind_gust_distribution_over_average_day(self):
         gs = fig.add_gridspec(nrows, ncols, hspace=0, wspace=0, left=0.12)
         axs = gs.subplots(sharex=True, sharey=True)
         x_min, x_max, y_min, y_max = 0, 0, 0, 0
-        for (i, j) in [(i, j) for i in range(nrows) for j in range(ncols)]:
+        for i, j in [(i, j) for i in range(nrows) for j in range(ncols)]:
             hour = hours[i * ncols + j]
             ax = axs[i, j]
             ax.set_title(f"{float(hour):05.2f}".replace(".", ":"), y=0.8)
@@ -2290,7 +2288,7 @@ def plot_sustained_surface_wind_speed_distribution_over_average_day(
         gs = fig.add_gridspec(nrows, ncols, hspace=0, wspace=0, left=0.12)
         axs = gs.subplots(sharex=True, sharey=True)
         x_min, x_max, y_min, y_max = 0, 0, 0, 0
-        for (i, j) in [(i, j) for i in range(nrows) for j in range(ncols)]:
+        for i, j in [(i, j) for i in range(nrows) for j in range(ncols)]:
             hour = hours[i * ncols + j]
             ax = axs[i, j]
             ax.set_title(f"{float(hour):05.2f}".replace(".", ":"), y=0.8)
@@ -2334,7 +2332,7 @@ def plot_sustained_surface_wind_speed_distribution_over_average_day(
         )
 
         if windSpeedLimit:
-            for (i, j) in [(i, j) for i in range(nrows) for j in range(ncols)]:
+            for i, j in [(i, j) for i in range(nrows) for j in range(ncols)]:
                 # Clear Sky Range Altitude Limits j]
                 ax = axs[i, j]
                 ax.vlines(
@@ -2656,7 +2654,7 @@ def plot_wind_profile_over_average_day(self, clear_range_limits=False):
         gs = fig.add_gridspec(nrows, ncols, hspace=0, wspace=0, left=0.12)
         axs = gs.subplots(sharex=True, sharey=True)
         x_min, x_max, y_min, y_max = 0, 0, np.inf, 0
-        for (i, j) in [(i, j) for i in range(nrows) for j in range(ncols)]:
+        for i, j in [(i, j) for i in range(nrows) for j in range(ncols)]:
             hour = hours[i * ncols + j]
             ax = axs[i, j]
             ax.plot(*self.average_wind_profile_at_given_hour[hour], "r-")
@@ -2689,7 +2687,7 @@ def plot_wind_profile_over_average_day(self, clear_range_limits=False):
         )
 
         if clear_range_limits:
-            for (i, j) in [(i, j) for i in range(nrows) for j in range(ncols)]:
+            for i, j in [(i, j) for i in range(nrows) for j in range(ncols)]:
                 # Clear Sky Range Altitude Limits
                 ax = axs[i, j]
                 ax.fill_between(
@@ -2794,7 +2792,7 @@ def plot_wind_heading_profile_over_average_day(self, clear_range_limits=False):
         gs = fig.add_gridspec(nrows, ncols, hspace=0, wspace=0, left=0.12)
         axs = gs.subplots(sharex=True, sharey=True)
         x_min, x_max, y_min, y_max = 0, 0, np.inf, 0
-        for (i, j) in [(i, j) for i in range(nrows) for j in range(ncols)]:
+        for i, j in [(i, j) for i in range(nrows) for j in range(ncols)]:
             hour = hours[i * ncols + j]
             ax = axs[i, j]
             ax.plot(*self.average_wind_heading_profile_at_given_hour[hour], "r-")
@@ -2819,7 +2817,7 @@ def plot_wind_heading_profile_over_average_day(self, clear_range_limits=False):
         )
 
         if clear_range_limits:
-            for (i, j) in [(i, j) for i in range(nrows) for j in range(ncols)]:
+            for i, j in [(i, j) for i in range(nrows) for j in range(ncols)]:
                 # Clear Sky range limits
                 ax = axs[i, j]
                 ax.fill_between(
diff --git a/rocketpy/Function.py b/rocketpy/Function.py
index f0d64e60b..4b9d09488 100644
--- a/rocketpy/Function.py
+++ b/rocketpy/Function.py
@@ -2085,7 +2085,7 @@ def __sub__(self, other):
                 # Create new Function object
                 return Function(source, inputs, outputs, interpolation)
             else:
-                return Function(lambda x: (self.getValue(x) * other(x)))
+                return Function(lambda x: (self.getValue(x) - other(x)))
         # If other is Float except...
         except AttributeError:
             if isinstance(other, (float, int, complex)):
diff --git a/rocketpy/__init__.py b/rocketpy/__init__.py
index bcc164608..c47529d39 100644
--- a/rocketpy/__init__.py
+++ b/rocketpy/__init__.py
@@ -16,7 +16,7 @@
 __copyright__ = "Copyright 20XX, Projeto Jupiter"
 __credits__ = ["Matheus Marques Araujo", "Rodrigo Schmitt", "Guilherme Tavares"]
 __license__ = "MIT"
-__version__ = "0.13.0"
+__version__ = "0.13.1"
 __maintainer__ = "Giovani Hidalgo Ceotto"
 __email__ = "ghceotto@gmail.com"
 __status__ = "Production"
diff --git a/setup.py b/setup.py
index 16feee477..591ae3447 100644
--- a/setup.py
+++ b/setup.py
@@ -1,16 +1,36 @@
+import os
+import sys
 import setuptools
 
 with open("README.md", "r") as fh:
     long_description = fh.read()
 
+# Manage different netCDF4 versions depending on system version and Google Colab
+netCDF4_requirement = "netCDF4>=1.6.2"
+if sys.version_info[1] == 7:
+    # Support for Python 3.7
+    netCDF4_requirement = "netCDF4>=1.4,<1.6"
+
+# Check if running on a linux system
+elif sys.platform == "linux":
+    if sys.version_info[1] < 11:
+        # Support for Linux with Python < 3.11
+        # Avoid problems due to netcdf4 issue linked below
+        netCDF4_requirement = "netCDF4>=1.4,<1.6"
+    else:
+        # Support for Linux with Python >= 3.11
+        # Might have problems with openDAP in the Environment class due to
+        # https://github.com/Unidata/netcdf4-python/issues/1179
+        netCDF4_requirement = "netCDF4>=1.6.2"
+
 setuptools.setup(
     name="rocketpy",
-    version="0.13.0",
+    version="0.13.1",
     install_requires=[
         "numpy>=1.0",
         "scipy>=1.0",
         "matplotlib>=3.0",
-        "netCDF4>=1.4,<1.6",
+        netCDF4_requirement,
         "windrose>=1.6.8",
         "ipywidgets>=7.6.3",
         "requests",
diff --git a/tests/fixtures/acceptance/EPFL_Bella_Lui/bella_lui_flight_sim.py b/tests/fixtures/acceptance/EPFL_Bella_Lui/bella_lui_flight_sim.py
index dfa7de298..cdfc52d45 100644
--- a/tests/fixtures/acceptance/EPFL_Bella_Lui/bella_lui_flight_sim.py
+++ b/tests/fixtures/acceptance/EPFL_Bella_Lui/bella_lui_flight_sim.py
@@ -112,6 +112,7 @@
     distanceToCM=parameters.get("tailDistanceToCM")[0],
 )
 
+
 # Parachute set-up
 def drogueTrigger(p, y):
     # p = pressure
diff --git a/tests/fixtures/acceptance/NDRT_2020/ndrt_2020_flight_sim.py b/tests/fixtures/acceptance/NDRT_2020/ndrt_2020_flight_sim.py
index 80b0f3373..854a0b981 100644
--- a/tests/fixtures/acceptance/NDRT_2020/ndrt_2020_flight_sim.py
+++ b/tests/fixtures/acceptance/NDRT_2020/ndrt_2020_flight_sim.py
@@ -122,6 +122,7 @@
     distanceToCM=parameters.get("transitiondistanceToCM")[0],
 )
 
+
 # Parachute set-up
 def drogueTrigger(p, y):
     # p = pressure