diff --git a/.readthedocs.yaml b/.readthedocs.yaml
new file mode 100644
index 0000000..be70604
--- /dev/null
+++ b/.readthedocs.yaml
@@ -0,0 +1,13 @@
+version: 2
+
+build:
+  os: ubuntu-22.04
+  tools:
+    python: "3.11"
+
+sphinx:
+  configuration: docs/conf.py
+
+python:
+  install:
+  - requirements: requirements.txt
\ No newline at end of file
diff --git a/docs/conf.py b/docs/conf.py
index 1c61bea..a3b84d3 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -13,7 +13,9 @@
     'sphinx.ext.autosummary',
     'sphinx.ext.mathjax',
     'numpydoc',
-    'sphinx.ext.graphviz'
+    'sphinx.ext.graphviz',
+	'sphinx.ext.viewcode',
+	'nbsphinx'
 ]
 
 source_suffix = '.rst'
@@ -29,6 +31,7 @@
 html_theme = 'nature'
 pygments_style = 'colorful'
 # htmlhelp_basename = 'JiTCODEdoc'
+exclude_patterns = ['_build', '**.ipynb_checkpoints']
 
 numpydoc_show_class_members = False
 autodoc_member_order = 'bysource'
diff --git a/docs/index.rst b/docs/index.rst
index 92e8e67..55e27d6 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -27,14 +27,16 @@ Indices and tables
 * :ref:`modindex`
 * :ref:`search`
 
-.. toctree::
-   :maxdepth: 2
-   :caption: Contents:
 
 
-Quick Start
+Notebooks
 *****************
 
+.. toctree::
+   :maxdepth: 2
+   :caption: Contents:
+   
+   notebooks
 
 Modules 
 *****************
diff --git a/docs/notebooks.rst b/docs/notebooks.rst
new file mode 100644
index 0000000..a1bc954
--- /dev/null
+++ b/docs/notebooks.rst
@@ -0,0 +1,9 @@
+
+Introduction
+............
+
+.. toctree::
+   :maxdepth: 2
+   :caption: Contents:
+
+mndynamics/examples/Intro
\ No newline at end of file
diff --git a/mndynamics/models/py/Erisir_Base.py b/mndynamics/models/py/Erisir_Base.py
index 1b184f5..fbce930 100644
--- a/mndynamics/models/py/Erisir_Base.py
+++ b/mndynamics/models/py/Erisir_Base.py
@@ -263,8 +263,7 @@ class Erisir_Burst(Erisir):
     The Erisir neuron with a slow potassium current that is strengthened by firing
     which turns the neuron into a burster.
 
-    Reference
-    ---------
+    Reference:
     An introduction to modeling neuronal dynamics, Borgers, Chapter 19.
     '''
 
diff --git a/mndynamics/models/py/FN_Base.py b/mndynamics/models/py/FN_Base.py
index db926cb..f088a97 100644
--- a/mndynamics/models/py/FN_Base.py
+++ b/mndynamics/models/py/FN_Base.py
@@ -60,17 +60,19 @@ def f_sys(self, x0, t):
 
     def simulate(self, tspan=None, x0=None):
         """
-        simulate the model
+        simulate the FN model
 
         Parameters
         ----------
-        tspan : array
-            time span for simulation
+        tspan : array-like, optional
+            time span for simulation, default is None
         
         Returns
         -------
         dict: {t, v, n}
-            time series of v, n
+            t: time span
+            v: membrane potential
+            n: gating variable
 
         """
         
diff --git a/mndynamics/models/py/HH_Base.py b/mndynamics/models/py/HH_Base.py
index 613fed8..3b7faaf 100644
--- a/mndynamics/models/py/HH_Base.py
+++ b/mndynamics/models/py/HH_Base.py
@@ -9,7 +9,6 @@ class HH(object):
     """
     Hudgkin Huxley Model
 
-    Usage:
     >>> par = {'i_ext': 1.5, 't_end': 100.0, 'v0': -70.0, 'dt': 0.01}
     >>> model = HH(par)
     >>> sol = model.simulate()
diff --git a/mndynamics/models/py/Inapik_Base.py b/mndynamics/models/py/Inapik_Base.py
index 5699f03..9a2dde6 100644
--- a/mndynamics/models/py/Inapik_Base.py
+++ b/mndynamics/models/py/Inapik_Base.py
@@ -11,7 +11,6 @@ class Inapik(object):
     """
     Inapik neuron model 
 
-    Usage:
     >>> par = {'}
     >>> model = Inapik(par)
     >>> sol = model.simulate()
diff --git a/mndynamics/models/py/LIF_Base.py b/mndynamics/models/py/LIF_Base.py
index b45fad4..031fb27 100644
--- a/mndynamics/models/py/LIF_Base.py
+++ b/mndynamics/models/py/LIF_Base.py
@@ -9,7 +9,6 @@ class LIF(object):
     """
     Linear Integrate-and-Fire (LIF) Model
 
-    Usage:
     >>> par = {'i_ext': 0.11, 't_end': 100.0, 'v0': -70.0, 'dt': 0.01}
     >>> model = LIF(par)
     >>> sol = model.simulate()
@@ -99,7 +98,6 @@ class LIF_Addapt(LIF):
     """
     Linear Integrate-and-Fire (LIF) Model with Adaptation
 
-    Usage:
     >>> par = {'i_ext': 0.13, 't_end': 100.0, 'v0': -70.0, 'dt': 0.01}
     >>> model = LIF_Addapt(par)
     >>> sol = model.simulate()
diff --git a/mndynamics/models/py/RTM_Base.py b/mndynamics/models/py/RTM_Base.py
index f3c0859..bf76b43 100644
--- a/mndynamics/models/py/RTM_Base.py
+++ b/mndynamics/models/py/RTM_Base.py
@@ -11,7 +11,6 @@ class RTM(object):
     """
     Reduced Traub-Miles Model of a Pyramidal Neuron in Rat Hippocampus
 
-    Usage:
     >>> par = {'i_ext': 1.5, 't_end': 100.0, 'v0': -70.0, 'dt': 0.01}
     >>> model = RTM(par)
     >>> sol = model.simulate()
@@ -132,7 +131,6 @@ class RTM_M(RTM):
     Reduced Traub-Miles neuron model with M-current
     Simulate the model with parameters given in the dictionary par:
 
-    Usage:
     >>> par = {'i_ext': 1.5, 't_end': 100.0, 'v0': -70.0, 'dt': 0.01}
     >>> model = RTM_M(par)
     >>> sol = model.simulate()
@@ -310,8 +308,7 @@ class RTM_2D(RTM):
     Reduced 2-dimensional Traub-Miles Model of a Pyramidal Neuron in Rat Hippocampus
     as defined in Eq. 12.1-2 in [1]
 
-    References
-    ----------
+    References:
     [1] Börgers, C., 2017. An introduction to modeling neuronal dynamics (Vol. 66). Berlin: Springer.
     '''
 
diff --git a/mndynamics/models/py/ch03.py b/mndynamics/models/py/ch03.py
index d6cc52c..7fbeed9 100644
--- a/mndynamics/models/py/ch03.py
+++ b/mndynamics/models/py/ch03.py
@@ -11,6 +11,10 @@ def __init__(self, par={}):
         super().__init__(par)
 
     def plot(self):
+        '''
+        Plot gating variables and time constants of the HH model
+        
+        '''
 
         fig, ax = plt.subplots(nrows=3, ncols=2, figsize=(7, 7))
         v = np.arange(-100, 50, 0.01)
diff --git a/mndynamics/models/py/ch04.py b/mndynamics/models/py/ch04.py
index 482fcec..1b6a32b 100644
--- a/mndynamics/models/py/ch04.py
+++ b/mndynamics/models/py/ch04.py
@@ -12,6 +12,10 @@ def __init__(self, par={}):
 
 
 class HH_REFRACTORINESS(_HH):
+    '''
+    Refractoriness of the HH model.
+    
+    '''
 
     def __init__(self, par=None):
         super().__init__(par)
@@ -33,6 +37,24 @@ def f_sys(self, x0, t, i_ext, PULSE_ONSET):
         return [dv, dm, dh, dn]
     
     def simulate(self, tspan=None, *args):
+        '''
+        Simulate the HH model.
+
+        Parameters
+        ----------
+        tspan : array_like, optional
+            Time points (ms) at which to solve for membrane potential and gating variables.
+            Default is np.arange(0, 50, 0.01).
+        args : tuple, optional
+            Additional arguments to pass to the ODE solver.
+
+        Returns
+        -------
+        dict
+            Dictionary with keys "t", "v", "m", "h", "n" and values being the corresponding
+            time points and membrane potential and gating variables.
+        
+        '''
         
         x0 = self.set_initial_state()