diff --git a/slycot/__init__.py b/slycot/__init__.py
index 69527783..32b4a5db 100644
--- a/slycot/__init__.py
+++ b/slycot/__init__.py
@@ -27,8 +27,9 @@
     # U : Utility Routines
 
 
-    # Analysis routines (16/60 wrapped)
+    # Analysis routines (17/60 wrapped)
     from .analysis import (ab01nd,
+                           ab04md,
                            ab05md, ab05nd,
                            ab07nd,
                            ab08nd, ab08nz,
diff --git a/slycot/analysis.py b/slycot/analysis.py
index 32050236..f84ff43f 100644
--- a/slycot/analysis.py
+++ b/slycot/analysis.py
@@ -139,6 +139,81 @@ def ab01nd(n, m, A, B, jobz='N', tol=0, ldwork=None):
         Z = None
     return Ac, Bc, ncont, indcon, nblk, Z, tau
 
+def ab04md(type_t, n, m, p, A, B, C, D, alpha=1.0, beta=1.0, ldwork=None):
+    """ At,Bt,Ct,Dt = ab04md(type_t, n, m, p, A, B, C, D, [alpha, beta,ldwork])
+
+    Parameters
+    ----------
+    type_t : {'D','C'}
+            Indicates the type of the original system and the
+            transformation to be performed as follows:
+            = 'D':  discrete-time   -> continuous-time;
+            = 'C':  continuous-time -> discrete-time.
+    n : int
+        The order of the matrix A, the number of rows of matrix B and
+        the number of columns of matrix C. It represents the dimension of
+        the state vector.  n > 0.        
+    m : int
+        The number of columns of matrix B. It represents the dimension of
+        the input vector.  m > 0.
+    p : int
+        The number of rows of matrix C. It represents the dimension of
+        the output vector.  p > 0.
+    A : (n,n) array_like
+        The leading n-by-n part of this array must contain the system state
+        matrix A.
+    B : (n,m) array_like
+        The leading n-by-m part of this array must contain the system input
+        matrix B.
+    C : (p,n) array_like
+        The leading p-by-n part of this array must contain the system output
+        matrix C.
+    D : (p,m) array_like
+        The leading p-by-m part of this array must contain the system direct
+        transmission matrix D.
+    alpha : double, optional
+        Parameter specifying the bilinear transformation.
+        Recommended values for stable systems: alpha = 1, alpha != 0,
+        Default is 1.0.
+    beta : double, optional
+        Parameter specifying the bilinear transformation.
+        Recommended values for stable systems: beta = 1, beta != 0,
+        Default is 1.0.
+    ldwork : int, optional
+        The length of the cache array.
+        ldwork >= max(1, n), default is max(1, n)
+    Returns
+    -------
+    At : (n,n) ndarray
+        The state matrix At of the transformed system.
+    Bt : (n,m) ndarray
+        The input matrix Bt of the transformed system.
+    Ct : (p,n) ndarray
+        The output matrix Ct of the transformed system.
+    Dt : (p,m) ndarray
+        The transmission matrix Dt of the transformed system.
+    Raises
+    ------
+    SlycotArithmeticError
+        :info == 1: 
+            If the matrix (ALPHA*I + A) is exactly singular
+        :info == 2: 
+            If the matrix  (BETA*I - A) is exactly singular.
+    """
+
+    hidden = ' (hidden by the wrapper)'  
+    arg_list = ['type_t', 'n', 'm', 'p', 'alpha', 'beta', 
+                'A', 'LDA'+hidden, 'B', 'LDB'+hidden, 'C', 'LDC'+hidden, 'D', 'LDD'+hidden,
+                'IWORK'+hidden, 'DWORK'+hidden, 'ldwork', 'info'+hidden]
+
+    if ldwork is None:
+        ldwork = max(1, n)
+
+    out = _wrapper.ab04md(type_t, n, m, p, alpha, beta, A, B, C, D, ldwork=ldwork)
+    info=out[-1]
+    raise_if_slycot_error(info, arg_list)
+
+    return out[:-1]
 
 def ab05md(n1,m1,p1,n2,p2,A1,B1,C1,D1,A2,B2,C2,D2,uplo='U'):
     """ n,a,b,c,d = ab05md(n1,m1,p1,n2,p2,a1,b1,c1,d1,a2,b2,c2,d2,[uplo])
diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf
index 633ff6ec..b2bb5a12 100644
--- a/slycot/src/analysis.pyf
+++ b/slycot/src/analysis.pyf
@@ -19,6 +19,26 @@ subroutine ab01nd(jobz,n,m,a,lda,b,ldb,ncont,indcon,nblk,z,ldz,tau,tol,iwork,dwo
     integer :: ldwork = max(n,3*m)
     integer intent(out) :: info
 end subroutine ab01nd
+subroutine ab04md(type_t,n,m,p,alpha,beta,a,lda,b,ldb,c,ldc,d,ldd,iwork,dwork,ldwork,info) ! in AB04MD.f
+    character :: type_t
+    integer check(n>=0) :: n
+    integer check(m>=0) :: m
+    integer check(p>=0) :: p
+    double precision intent(in) :: alpha
+    double precision intent(in) :: beta
+    double precision intent(in,out,copy), dimension(n,n),depend(n) :: a
+    integer intent(hide),depend(a) :: lda = shape(a,0)
+    double precision intent(in,out,copy), dimension(n,m),depend(n,m) :: b
+    integer intent(hide),depend(b) :: ldb = shape(b,0)
+    double precision intent(in,out,copy), dimension(p,n),depend(n,p) :: c
+    integer intent(hide),depend(c) :: ldc = shape(c,0)
+    double precision intent(in,out,copy), dimension(p,m),depend(m,p) :: d
+    integer intent(hide),depend(d) :: ldd = shape(d,0)
+    integer intent(hide,cache),dimension(n),depend(n) :: iwork
+    double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
+    integer :: ldwork = max(1,n)
+    integer intent(out) :: info
+end subroutine ab04md
 subroutine ab05md(uplo,over,n1,m1,p1,n2,p2,a1,lda1,b1,ldb1,c1,ldc1,d1,ldd1,a2,lda2,b2,ldb2,c2,ldc2,d2,ldd2,n,a,lda,b,ldb,c,ldc,d,ldd,dwork,ldwork,info) ! in AB05MD.f
     character :: uplo = 'U'
     character intent(hide) :: over = 'N' ! not sure how the overlap works
diff --git a/slycot/tests/CMakeLists.txt b/slycot/tests/CMakeLists.txt
index 2aafab3f..80d751da 100644
--- a/slycot/tests/CMakeLists.txt
+++ b/slycot/tests/CMakeLists.txt
@@ -2,6 +2,7 @@ set(PYSOURCE
 
   __init__.py
   test_ab01.py
+  test_ab04md.py
   test_ab08n.py
   test_ag08bd.py
   test_examples.py
diff --git a/slycot/tests/test_ab04md.py b/slycot/tests/test_ab04md.py
new file mode 100644
index 00000000..7ce59a0e
--- /dev/null
+++ b/slycot/tests/test_ab04md.py
@@ -0,0 +1,66 @@
+from slycot import analysis
+import numpy as np
+
+from numpy.testing import assert_allclose
+
+
+class test_ab04md:
+    """Test ab04md.
+
+    Example data taken from
+    https://www.slicot.org/objects/software/shared/doc/AB04MD.html
+    """
+
+    Ac = np.array([[1.0, 0.5],
+                   [0.5, 1.0]])
+    Bc = np.array([[0.0, -1.0],
+                   [1.0, 0.0]])
+    Cc = np.array([[-1.0, 0.0],
+                   [0.0, 1.0]])
+    Dc = np.array([[1.0, 0.0],
+                   [0.0, -1.0]])
+
+    Ad = np.array([[-1.0, -4.0],
+                   [-4.0, -1.0]])
+    Bd = np.array([[2.8284, 0.0],
+                   [0.0, -2.8284]])
+    Cd = np.array([[0.0, 2.8284],
+                   [-2.8284, 0.0]])
+    Dd = np.array([[-1.0, 0.0],
+                   [0.0, -3.0]])
+
+    def test_ab04md_cont_disc_cont(self):
+        """Test transformation from continuous - to discrete - to continuous time.
+        """
+
+        n, m = self.Bc.shape
+        p = self.Cc.shape[0]
+
+        Ad_t, Bd_t, Cd_t, Dd_t = analysis.ab04md(
+            'C', n, m, p, self.Ac, self.Bc, self.Cc, self.Dc)
+
+        Ac_t, Bc_t, Cc_t, Dc_t = analysis.ab04md(
+            'D', n, m, p, Ad_t, Bd_t, Cd_t, Dd_t)
+
+        assert_allclose(self.Ac, Ac_t)
+        assert_allclose(self.Bc, Bc_t)
+        assert_allclose(self.Cc, Cc_t)
+        assert_allclose(self.Dc, Dc_t)
+
+    def test_ab04md_disc_cont_disc(self):
+        """Test transformation from discrete - to continuous - to discrete time.
+        """
+
+        n, m = self.Bc.shape
+        p = self.Cc.shape[0]
+
+        Ac_t, Bc_t, Cc_t, Dc_t = analysis.ab04md(
+            'D', n, m, p, self.Ad, self.Bd, self.Cd, self.Dd)
+
+        Ad_t, Bd_t, Cd_t, Dd_t = analysis.ab04md(
+            'C', n, m, p, Ac_t, Bc_t, Cc_t, Dc_t)
+
+        assert_allclose(self.Ad, Ad_t)
+        assert_allclose(self.Bd, Bd_t)
+        assert_allclose(self.Cd, Cd_t)
+        assert_allclose(self.Dd, Dd_t)