diff --git a/CHANGES b/CHANGES
--- a/CHANGES
+++ b/CHANGES
@@ -1,3 +1,12 @@
+0.8.3.0
+=======
+
+- odeSolve
+
+- Matrix arithmetic automatically replicates matrix with single row/column
+
+- latexFormat, dispcf
+
 0.8.2.0
 =======
 
@@ -56,3 +65,4 @@
 - added NFData instances for Matrix and Vector.
 
 - liftVector, liftVector2 replaced by mapVector, zipVector.
+
diff --git a/README b/README
--- a/README
+++ b/README
@@ -50,8 +50,8 @@
 
 ACKNOWLEDGEMENTS -----------------------------------------------------
 
-I thank Don Stewart, Henning Thielemann, Bulat Ziganshin and all the people
-in the Haskell mailing lists for their help.
+I thank Don Stewart, Henning Thielemann, Bulat Ziganshin, Heinrich Apfelmus,
+and all the people in the Haskell mailing lists for their help.
 
 - Nico Mahlo discovered a bug in the eigendecomposition wrapper.
 
@@ -108,3 +108,5 @@
   to the values at that position.
 
 - Jean-Francois Tremblay discovered an error in the tutorial.
+
+- Heinrich Apfelmus fixed hmatrix.cabal for OS/X.
diff --git a/examples/Real.hs b/examples/Real.hs
new file mode 100644
--- /dev/null
+++ b/examples/Real.hs
@@ -0,0 +1,104 @@
+
+-- Alternative interface and utilities for creation of real arrays, useful to work in interactive mode.
+
+module Real(
+    module Numeric.LinearAlgebra,
+    (<>), (*>), (<*), (<\>), (\>),
+    vector,
+    eye,
+    zeros, ones,
+    diagl,
+    row,
+    col,
+    (#),(&), (//), blocks,
+    rand
+) where
+
+import Numeric.LinearAlgebra hiding ((<>), (<|>), (<->), (<\>), (.*), (*/))
+import System.Random(randomIO)
+
+infixl 7 <>
+-- | Matrix product ('multiply')
+(<>) :: Field t => Matrix t -> Matrix t -> Matrix t
+(<>) = multiply
+
+infixl 7 *>
+-- | matrix x vector
+(*>) :: Field t => Matrix t -> Vector t -> Vector t
+m *> v = flatten $ m <> (asColumn v)
+
+infixl 7 <*
+-- | vector x matrix
+(<*) :: Field t => Vector t -> Matrix t -> Vector t
+v <* m = flatten $ (asRow v) <> m
+
+
+-- | Least squares solution of a linear system for several right-hand sides, similar to the \\ operator of Matlab\/Octave. (\<\\\>) = 'linearSolveSVD'.
+(<\>) :: (Field a) => Matrix a -> Matrix a -> Matrix a
+infixl 7 <\>
+(<\>) = linearSolveSVD
+
+-- | Least squares solution of a linear system for a single right-hand side. See '(\<\\\>)'.
+(\>) :: (Field a) => Matrix a -> Vector a -> Vector a
+infixl 7 \>
+m \> v = flatten (m <\> reshape 1 v)
+
+-- | Pseudorandom matrix with uniform elements between 0 and 1.
+rand :: Int -- ^ rows
+     -> Int -- ^ columns
+     -> IO (Matrix Double)
+rand r c = do
+    seed <- randomIO
+    return (reshape c $ randomVector seed Uniform (r*c))
+
+-- | Real identity matrix.
+eye :: Int -> Matrix Double
+eye = ident
+
+-- | Create a real vector from a list.
+vector :: [Double] -> Vector Double
+vector = fromList
+
+-- | Create a real diagonal matrix from a list.
+diagl :: [Double] -> Matrix Double
+diagl = diag . vector
+
+-- | Create a matrix or zeros.
+zeros :: Int -- ^ rows
+      -> Int -- ^ columns
+      -> Matrix Double
+zeros r c = reshape c (constant 0 (r*c))
+
+-- | Create a matrix or ones.
+ones :: Int -- ^ rows
+     -> Int -- ^ columns
+     -> Matrix Double
+ones r c = reshape c (constant 1 (r*c))
+
+-- | Concatenation of real vectors.
+infixl 9 #
+(#) :: Vector Double -> Vector Double -> Vector Double
+a # b = join [a,b]
+
+-- | Horizontal concatenation of real matrices.
+infixl 8 &
+(&) :: Matrix Double -> Matrix Double -> Matrix Double
+a & b = fromBlocks [[a,b]]
+
+-- | Vertical concatenation of real matrices.
+infixl 7 //
+(//) :: Matrix Double -> Matrix Double -> Matrix Double
+a // b = fromBlocks [[a],[b]]
+
+-- | Real block matrix from a rectangular list of lists.
+blocks :: [[Matrix Double]] -> Matrix Double
+blocks = fromBlocks
+
+-- | A real matrix with a single row, create from a list of elements.
+row :: [Double] -> Matrix Double
+row = asRow . vector
+
+-- | A real matrix with a single column, created from a list of elements.
+col :: [Double] -> Matrix Double
+col = asColumn . vector
+
diff --git a/examples/benchmarks.hs b/examples/benchmarks.hs
deleted file mode 100644
--- a/examples/benchmarks.hs
+++ /dev/null
@@ -1,151 +0,0 @@
-{-# LANGUAGE BangPatterns #-}
-
--- $ ghc --make -O2 benchmarks.hs
-
-
-import Numeric.LinearAlgebra
-import System.Time
-import System.CPUTime
-import Text.Printf
-import Data.List(foldl1')
-
-
-time act = do
-    t0 <- getCPUTime
-    act
-    t1 <- getCPUTime
-    printf "%.3f s CPU\n" $ (fromIntegral (t1 - t0) / (10^12 :: Double)) :: IO ()
-
---------------------------------------------------------------------------------
-
-main = sequence_ [bench1,
-                  bench2,
-                  bench4,
-                  bench5 1000000 3, bench5 100000 50,
-                  bench6 100 (100000::Double), bench6 100000 (100::Double), bench6 10000 (1000::Double)]
-
-w :: Vector Double
-w = constant 1 5000000
-w2 = 1 * w
-
-v = flatten $ ident 500 :: Vector Double
-
-
-bench1 = do
-    time $ print$ vectorMax (w+w2) -- evaluate it
-    putStrLn "Sum of a vector with 5M doubles:"
-    print$ vectorMax (w+w2) -- evaluate it
-    time $ printf "         BLAS: %.2f: " $ sumVB w
-    time $ printf "BLAS only dot: %.2f: " $ w <.> w2
-    time $ printf "      Haskell: %.2f: " $ sumVH w
-    time $ printf "       innerH: %.2f: " $ innerH w w2
-    time $ printf "foldVector: %.2f: " $ sumVector w
-    let getPos k s = if k `mod` 500 < 200 && w@>k > 0 then k:s else s
-    putStrLn "foldLoop for element selection:"
-    time $ print $ (`divMod` 500) $ maximum $ foldLoop getPos [] (dim w)
-    putStrLn "constant 5M:"
-    time $ print $ constant (1::Double) 5000001 @> 7
-    time $ print $ constant           i 5000001 @> 7
-    time $ print $ conj (constant i 5000001) @> 7
-    putStrLn "zips C vs H:"
-    time $ print $ (w / w2) @> 7
-    time $ print $ (zipVector (/) w w2) @> 7
-    putStrLn "folds C/BLAS vs H:"
-    let t = constant (1::Double) 5000002
-    print $ t @> 7
-    time $ print $ foldVector max (t@>0) t
-    time $ print $ vectorMax t
-    time $ print $ sqrt $ foldVector (\v s -> v*v+s) 0 t
-    time $ print $ pnorm PNorm2 t
-    putStrLn "scale C/BLAS vs H:"
-    time $ print $ mapVector (*2) t @> 7
-    time $ print $ (2 * t) @> 7
-
-sumVB v = constant 1 (dim v) <.> v
-
-sumVH v = go (d - 1) 0
-     where
-       d = dim v
-       go :: Int -> Double -> Double
-       go 0 s = s + (v @> 0)
-       go !j !s = go (j - 1) (s + (v @> j))
-
-innerH u v = go (d - 1) 0
-     where
-       d = min (dim u) (dim v)
-       go :: Int -> Double -> Double
-       go 0 s = s + (u @> 0) * (v @> 0)
-       go !j !s = go (j - 1) (s + (u @> j) * (v @> j))
-
-
--- sumVector = foldVectorG (\k v s -> v k + s) 0.0
-sumVector = foldVector (+) 0.0
-
---------------------------------------------------------------------------------
-
-bench2 = do
-    putStrLn "-------------------------------------------------------"
-    putStrLn "Multiplication of 1M different 3x3 matrices:"
---    putStrLn "from [[]]"
---    time $ print $ manymult (10^6) rot'
---    putStrLn "from (3><3) []"
-    time $ print $ manymult (10^6) rot
-    print $ cos (10^6/2)
-
-
-rot' :: Double -> Matrix Double
-rot' a = matrix [[ c,0,s],
-                 [ 0,1,0],
-                 [-s,0,c]]
-    where c = cos a
-          s = sin a
-          matrix = fromLists
-
-rot :: Double -> Matrix Double
-rot a = (3><3) [ c,0,s
-               , 0,1,0
-               ,-s,0,c ]
-    where c = cos a
-          s = sin a
-
-manymult n r = foldl1' (<>) (map r angles)
-    where angles = toList $ linspace n (0,1)
-          -- angles = map (k*) [0..n']
-          -- n' = fromIntegral n - 1
-          -- k  = recip n'
-
---------------------------------------------------------------------------------
-
-bench4 = do
-    putStrLn "-------------------------------------------------------"
-    putStrLn "1000x1000 inverse"
-    let a = ident 1000 :: Matrix Double
-    let b = 2*a
-    print $ vectorMax $ flatten (a+b) -- evaluate it
-    time $ print $ vectorMax $ flatten $ linearSolve a b
-
---------------------------------------------------------------------------------
-
-op1 a b = a <> trans b
-op2 a b = a  + trans b
-
-timep = time . print . vectorMax . flatten
-
-bench5 n d = do
-    putStrLn "-------------------------------------------------------"
-    putStrLn "transpose in add"
-    let ms = replicate n ((ident d :: Matrix Double))
-    timep $ foldl1' (+)  ms
-    timep $ foldl1' op2  ms
-    putStrLn "-------------------------------------------------------"
-    putStrLn "transpose in multiply"
-
-    timep $ foldl1' (<>) ms
-    timep $ foldl1' op1  ms
-
---------------------------------------------------------------------------------
-
-bench6 sz n = do
-    putStrLn "-------------------------------------------------------"
-    putStrLn "many constants"
-    time $ print $ sum $ map ((@>0). flip constant sz) [1..n]
diff --git a/examples/latexmat.hs b/examples/latexmat.hs
deleted file mode 100644
--- a/examples/latexmat.hs
+++ /dev/null
@@ -1,11 +0,0 @@
-import Numeric.LinearAlgebra
-import Text.Printf
-
-disp w l fmt m = unlines $ map (++l) $ lines $ format w (printf fmt) m
-
-latex fmt m = "\\begin{bmatrix}\n" ++ disp " & " " \\\\" fmt m ++ "\\end{bmatrix}"
-
-main = do
-    let m = (3><4) [1..12::Double]
-    putStrLn $ disp " | " "" "%.2f" m
-    putStrLn $ latex "%.3f" m
diff --git a/examples/ode.hs b/examples/ode.hs
new file mode 100644
--- /dev/null
+++ b/examples/ode.hs
@@ -0,0 +1,34 @@
+import Numeric.GSL.ODE
+import Numeric.LinearAlgebra
+import Graphics.Plot
+
+vanderpol mu = do
+    let xdot mu t [x,v] = [v, -x + mu * v * (1-x^2)]
+        ts = linspace 1000 (0,50)
+        sol = toColumns $ odeSolve (xdot mu) [1,0] ts
+    mplot (ts : sol)
+    mplot sol
+
+
+harmonic w d = do
+    let xdot w d t [x,v] = [v, a*x + b*v] where a = -w^2; b = -2*d*w
+        ts = linspace 100 (0,20)
+        sol = odeSolve (xdot w d) [1,0] ts
+    mplot (ts : toColumns sol)
+
+
+kepler v a = mplot (take 2 $ toColumns sol) where
+    xdot t [x,y,vx,vy] = [vx,vy,x*k,y*k]
+        where g=1
+              k=(-g)*(x*x+y*y)**(-1.5)
+    ts = linspace 100 (0,30)
+    sol = odeSolve xdot [4, 0, v * cos (a*degree), v * sin (a*degree)] ts
+    degree = pi/180
+
+
+main = do
+    vanderpol 2
+    harmonic 1 0
+    harmonic 1 0.1
+    kepler 0.3 60
+    kepler 0.4 70
diff --git a/examples/pca1.hs b/examples/pca1.hs
--- a/examples/pca1.hs
+++ b/examples/pca1.hs
@@ -8,16 +8,15 @@
 type Vec = Vector Double
 type Mat = Matrix Double
 
-sumColumns m = constant 1 (rows m) <> m
 
--- Vec with the mean value of the columns of a Mat
-mean x = sumColumns x / fromIntegral (rows x)
+-- Vector with the mean value of the columns of a matrix
+mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a
 
--- covariance Mat of a list of observations as rows of a Mat
-cov x = (trans xc <> xc) / fromIntegral (rows x -1)
-    where xc = center x
-          center m = m - constant 1 (rows m) `outer` mean m
+-- covariance matrix of a list of observations stored as rows
+cov x = (trans xc <> xc) / fromIntegral (rows x - 1)
+    where xc = x - asRow (mean x)
 
+
 -- creates the compression and decompression functions from the desired number of components
 pca :: Int -> Mat -> (Vec -> Vec , Vec -> Vec)
 pca n dataSet = (encode,decode)
@@ -38,7 +37,7 @@
         system("wget -nv http://dis.um.es/~alberto/material/sp/mnist.txt.gz")
         system("gunzip mnist.txt.gz")
         return ()
-    m <- fromFile "mnist.txt" (5000,785)
+    m <- loadMatrix "mnist.txt" -- fromFile "mnist.txt" (5000,785)
     let xs = takeColumns (cols m -1) m -- the last column is the digit type (class label)
     let x = toRows xs !! 4  -- an arbitrary test Vec
     let (pe,pd) = pca 10 xs
diff --git a/examples/pca2.hs b/examples/pca2.hs
--- a/examples/pca2.hs
+++ b/examples/pca2.hs
@@ -9,33 +9,31 @@
 type Vec = Vector Double
 type Mat = Matrix Double
 
-sumColumns m = constant 1 (rows m) <> m
+-- Vector with the mean value of the columns of a matrix
+mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a
 
--- Vector with the mean value of the columns of a Mat
-mean x = sumColumns x / fromIntegral (rows x)
+-- covariance matrix of a list of observations stored as rows
+cov x = (trans xc <> xc) / fromIntegral (rows x - 1)
+    where xc = x - asRow (mean x)
 
--- covariance matrix of a list of observations as rows of a matrix
-cov x = (trans xc <> xc) / fromIntegral (rows x -1) 
-    where xc = center x
-          center m = m - constant 1 (rows m) `outer` mean m
 
 type Stat = (Vec, [Double], Mat)
--- 1st and 2nd order statistics of a dataset (mean, eigenvalues and eigenvectors of cov)  
+-- 1st and 2nd order statistics of a dataset (mean, eigenvalues and eigenvectors of cov)
 stat :: Mat -> Stat
-stat x = (m, toList s, trans v) where   
+stat x = (m, toList s, trans v) where
     m = mean x
-    (s,v) = eigSH' (cov x)   
+    (s,v) = eigSH' (cov x)
 
--- creates the compression and decompression functions from the desired reconstruction 
+-- creates the compression and decompression functions from the desired reconstruction
 -- quality and the statistics of a data set
 pca :: Double -> Stat -> (Vec -> Vec , Vec -> Vec)
 pca prec (m,s,v) = (encode,decode)    
   where    
     encode x = vp <> (x - m)
     decode x = x <> vp + m
-    vp = takeRows n v    
+    vp = takeRows n v
     n = 1 + (length $ fst $ span (< (prec'*sum s)) $ cumSum s)
-    cumSum = tail . scanl (+) 0.0     
+    cumSum = tail . scanl (+) 0.0
     prec' = if prec <=0.0 || prec >= 1.0
                 then error "the precision in pca must be 0<prec<1"
                 else prec
@@ -49,7 +47,7 @@
     let (pe,pd) = pca prec st
     let y = pe x
     print $ dim y
-    shdigit (pd y)    
+    shdigit (pd y)
 
 main = do
     ok <- doesFileExist ("mnist.txt")
@@ -58,7 +56,7 @@
         system("wget -nv http://dis.um.es/~alberto/material/sp/mnist.txt.gz")
         system("gunzip mnist.txt.gz")
         return ()
-    m <- fromFile "mnist.txt" (5000,785)
+    m <- loadMatrix "mnist.txt"
     let xs = takeColumns (cols m -1) m
     let x = toRows xs !! 4  -- an arbitrary test vector
     shdigit x
diff --git a/examples/pinv.hs b/examples/pinv.hs
--- a/examples/pinv.hs
+++ b/examples/pinv.hs
@@ -3,7 +3,7 @@
 import Text.Printf(printf)
 
 expand :: Int -> Vector Double -> Matrix Double
-expand n x = fromColumns $ constant 1 (dim x): map (x^) [1 .. n]
+expand n x = fromColumns $ map (x^) [0 .. n]
 
 polynomialModel :: Vector Double -> Vector Double -> Int
                 -> (Vector Double -> Vector Double)
diff --git a/hmatrix.cabal b/hmatrix.cabal
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix
-Version:            0.8.2.0
+Version:            0.8.3.1
 License:            GPL
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -14,10 +14,11 @@
 tested-with:        GHC ==6.10.4, GHC ==6.12.1
 
 cabal-version:      >=1.2
+
 build-type:         Custom
-extra-source-files: lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h
 
-extra-source-files: configure configure.hs README INSTALL CHANGES
+extra-source-files: lib/Numeric/LinearAlgebra/Tests/quickCheckCompat.h
+                    configure configure.hs README INSTALL CHANGES
 extra-tmp-files:    hmatrix.buildinfo
 
 extra-source-files: examples/tests.hs
@@ -25,6 +26,7 @@
                     examples/integrate.hs
                     examples/minimize.hs
                     examples/root.hs
+                    examples/ode.hs
                     examples/pca1.hs
                     examples/pca2.hs
                     examples/pinv.hs
@@ -33,12 +35,11 @@
                     examples/kalman.hs
                     examples/parallel.hs
                     examples/plot.hs
-                    examples/latexmat.hs
                     examples/inplace.hs
-                    examples/benchmarks.hs
                     examples/error.hs
                     examples/devel/wrappers.hs
                     examples/devel/functions.c
+                    examples/Real.hs
 
 extra-source-files: lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h,
                     lib/Numeric/LinearAlgebra/LAPACK/clapack.h
@@ -47,28 +48,19 @@
                     lib/Numeric/GSL/Special/autoall.sh,
                     lib/Numeric/GSL/Special/replace.hs
 
-flag splitBase
-    description:    Choose the new smaller, split-up base package.
-
 flag mkl
     description:    Link with Intel's MKL optimized libraries.
     default:        False
 
-flag accelerate
-    description:    Use the accelerate framework for LAPACK/BLAS on OS X
-    default:        False
-
 flag unsafe
     description:    Compile the library with bound checking disabled.
     default:        False
 
 library
-    if flag(splitBase)
-      build-depends:    base >= 3 && < 5, array
-    else
-      build-depends:    base < 3
 
-    Build-Depends:      haskell98,
+    Build-Depends:      base >= 3 && < 5,
+                        array,
+                        haskell98,
                         QuickCheck, HUnit,
                         storable-complex,
                         process
@@ -86,6 +78,7 @@
                         Numeric.GSL.Polynomials,
                         Numeric.GSL.Minimization,
                         Numeric.GSL.Root,
+                        Numeric.GSL.ODE,
                         Numeric.GSL.Vector,
                         Numeric.GSL.Special,
                         Numeric.GSL.Special.Gamma,
@@ -127,7 +120,7 @@
                         Numeric.LinearAlgebra.Tests,
                         Data.Packed.Convert,
                         Data.Packed.ST,
-                        Data.Packed.Development
+                        Data.Packed.Development,
                         Data.Packed.Random
     other-modules:      Data.Packed.Internal,
                         Data.Packed.Internal.Common,
@@ -159,12 +152,19 @@
       else
         extra-libraries:   gsl mkl_lapack mkl_intel mkl_sequential mkl_core
 
-    if flag(accelerate)
-        frameworks: Accelerate
+    if os(OSX)
+        extra-lib-dirs: /opt/local/lib/
+        include-dirs: /opt/local/include
         extra-libraries: gsl
+        frameworks: Accelerate
 
 -- The extra-libraries required for GSL and LAPACK
 -- should now be automatically detected by configure(.hs)
 
     extra-libraries:
     extra-lib-dirs:
+
+    source-repository head
+        type:     darcs
+        location: http://code.haskell.org/hmatrix
+
diff --git a/lib/Data/Packed/Internal/Common.hs b/lib/Data/Packed/Internal/Common.hs
--- a/lib/Data/Packed/Internal/Common.hs
+++ b/lib/Data/Packed/Internal/Common.hs
@@ -18,8 +18,9 @@
   Adapt,
   app1, app2, app3, app4,
   (//), check,
-  partit, common, compatdim,
-  fi
+  splitEvery, common, compatdim,
+  fi,
+  table
 ) where
 
 import Foreign
@@ -27,11 +28,12 @@
 import Foreign.C.String(peekCString)
 import Foreign.C.Types
 import Foreign.Storable.Complex()
+import Data.List(transpose,intersperse)
 
--- | @partit 3 [1..9] == [[1,2,3],[4,5,6],[7,8,9]]@
-partit :: Int -> [a] -> [[a]]
-partit _ [] = []
-partit n l  = take n l : partit n (drop n l)
+-- | @splitEvery 3 [1..9] == [[1,2,3],[4,5,6],[7,8,9]]@
+splitEvery :: Int -> [a] -> [[a]]
+splitEvery _ [] = []
+splitEvery k l = take k l : splitEvery k (drop k l)
 
 -- | obtains the common value of a property of a list
 common :: (Eq a) => (b->a) -> [b] -> Maybe a
@@ -46,6 +48,15 @@
 compatdim [] = Nothing
 compatdim [a] = Just a
 compatdim (a:b:xs) = if a==b || a==1 || b==1 then compatdim (max a b:xs) else Nothing
+
+-- | Formatting tool
+table :: String -> [[String]] -> String
+table sep as = unlines . map unwords' $ transpose mtp where 
+    mt = transpose as
+    longs = map (maximum . map length) mt
+    mtp = zipWith (\a b -> map (pad a) b) longs mt
+    pad n str = replicate (n - length str) ' ' ++ str
+    unwords' = concat . intersperse sep
 
 -- | postfix function application (@flip ($)@)
 (//) :: x -> (x -> y) -> y
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -16,7 +16,7 @@
 -- #hide
 
 module Data.Packed.Internal.Matrix(
-    Matrix(..),
+    Matrix(..), rows, cols,
     MatrixOrder(..), orderOf,
     createMatrix, withMatrix, mat,
     cmat, fmat,
@@ -77,17 +77,23 @@
 data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq)
 
 -- | Matrix representation suitable for GSL and LAPACK computations.
-data Matrix t = MC { rows :: {-# UNPACK #-} !Int
-                   , cols :: {-# UNPACK #-} !Int
+data Matrix t = MC { irows :: {-# UNPACK #-} !Int
+                   , icols :: {-# UNPACK #-} !Int
                    , cdat :: {-# UNPACK #-} !(Vector t) }
 
-              | MF { rows :: {-# UNPACK #-} !Int
-                   , cols :: {-# UNPACK #-} !Int
+              | MF { irows :: {-# UNPACK #-} !Int
+                   , icols :: {-# UNPACK #-} !Int
                    , fdat :: {-# UNPACK #-} !(Vector t) }
 
 -- MC: preferred by C, fdat may require a transposition
 -- MF: preferred by LAPACK, cdat may require a transposition
 
+rows :: Matrix t -> Int
+rows = irows
+
+cols :: Matrix t -> Int
+cols = icols
+
 xdat MC {cdat = d } = d
 xdat MF {fdat = d } = d
 
@@ -97,16 +103,16 @@
 
 -- | Matrix transpose.
 trans :: Matrix t -> Matrix t
-trans MC {rows = r, cols = c, cdat = d } = MF {rows = c, cols = r, fdat = d }
-trans MF {rows = r, cols = c, fdat = d } = MC {rows = c, cols = r, cdat = d }
+trans MC {irows = r, icols = c, cdat = d } = MF {irows = c, icols = r, fdat = d }
+trans MF {irows = r, icols = c, fdat = d } = MC {irows = c, icols = r, cdat = d }
 
 cmat :: (Element t) => Matrix t -> Matrix t
 cmat m@MC{} = m
-cmat MF {rows = r, cols = c, fdat = d } = MC {rows = r, cols = c, cdat = transdata r d c}
+cmat MF {irows = r, icols = c, fdat = d } = MC {irows = r, icols = c, cdat = transdata r d c}
 
 fmat :: (Element t) => Matrix t -> Matrix t
 fmat m@MF{} = m
-fmat MC {rows = r, cols = c, cdat = d } = MF {rows = r, cols = c, fdat = transdata c d r}
+fmat MC {irows = r, icols = c, cdat = d } = MF {irows = r, icols = c, fdat = transdata c d r}
 
 -- C-Haskell matrix adapter
 mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r
@@ -133,7 +139,7 @@
 
 -- | the inverse of 'Data.Packed.Matrix.fromLists'
 toLists :: (Element t) => Matrix t -> [[t]]
-toLists m = partit (cols m) . toList . flatten $ m
+toLists m = splitEvery (cols m) . toList . flatten $ m
 
 -- | Create a matrix from a list of vectors.
 -- All vectors must have the same dimension,
@@ -170,13 +176,13 @@
 --    | i<0 || i>=r || j<0 || j>=c = error "matrix indexing out of range"
 --    | otherwise   = cdat m `at` (i*c+j)
 
-MC {rows = r, cols = c, cdat = v} @@> (i,j)
+MC {irows = r, icols = c, cdat = v} @@> (i,j)
     | safe      = if i<0 || i>=r || j<0 || j>=c
                     then error "matrix indexing out of range"
                     else v `at` (i*c+j)
     | otherwise = v `at` (i*c+j)
 
-MF {rows = r, cols = c, fdat = v} @@> (i,j)
+MF {irows = r, icols = c, fdat = v} @@> (i,j)
     | safe      = if i<0 || i>=r || j<0 || j>=c
                     then error "matrix indexing out of range"
                     else v `at` (j*r+i)
@@ -184,18 +190,18 @@
 {-# INLINE (@@>) #-}
 
 --  Unsafe matrix access without range checking
-atM' MC {cols = c, cdat = v} i j = v `at'` (i*c+j)
-atM' MF {rows = r, fdat = v} i j = v `at'` (j*r+i)
+atM' MC {icols = c, cdat = v} i j = v `at'` (i*c+j)
+atM' MF {irows = r, fdat = v} i j = v `at'` (j*r+i)
 {-# INLINE atM' #-}
 
 ------------------------------------------------------------------
 
-matrixFromVector RowMajor c v = MC { rows = r, cols = c, cdat = v }
+matrixFromVector RowMajor c v = MC { irows = r, icols = c, cdat = v }
     where (d,m) = dim v `divMod` c
           r | m==0 = d
             | otherwise = error "matrixFromVector"
 
-matrixFromVector ColumnMajor c v = MF { rows = r, cols = c, fdat = v }
+matrixFromVector ColumnMajor c v = MF { irows = r, icols = c, fdat = v }
     where (d,m) = dim v `divMod` c
           r | m==0 = d
             | otherwise = error "matrixFromVector"
@@ -223,8 +229,8 @@
 
 -- | application of a vector function on the flattened matrix elements
 liftMatrix :: (Element a, Element b) => (Vector a -> Vector b) -> Matrix a -> Matrix b
-liftMatrix f MC { cols = c, cdat = d } = matrixFromVector RowMajor    c (f d)
-liftMatrix f MF { cols = c, fdat = d } = matrixFromVector ColumnMajor c (f d)
+liftMatrix f MC { icols = c, cdat = d } = matrixFromVector RowMajor    c (f d)
+liftMatrix f MF { icols = c, fdat = d } = matrixFromVector ColumnMajor c (f d)
 
 -- | application of a vector function on the flattened matrices elements
 liftMatrix2 :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -15,7 +15,7 @@
 -- #hide
 
 module Data.Packed.Internal.Vector (
-    Vector(..),
+    Vector(..), dim,
     fromList, toList, (|>),
     join, (@>), safe, at, at', subVector,
     mapVector, zipVector,
@@ -47,10 +47,14 @@
 
 -- | A one-dimensional array of objects stored in a contiguous memory block.
 data Vector t =
-    V { dim  :: {-# UNPACK #-} !Int               -- ^ number of elements
+    V { idim  :: {-# UNPACK #-} !Int              -- ^ number of elements
       , fptr :: {-# UNPACK #-} !(ForeignPtr t)    -- ^ foreign pointer to the memory block
       }
 
+-- | Number of elements
+dim :: Vector t -> Int
+dim = idim
+
 -- C-Haskell vector adapter
 vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r
 vec = withVector
@@ -157,7 +161,7 @@
                         -> Int       -- ^ number of elements to extract
                         -> Vector t  -- ^ source
                         -> Vector t  -- ^ result
-subVector k l (v@V {dim=n})
+subVector k l (v@V {idim=n})
     | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range"
     | otherwise = unsafePerformIO $ do
         r <- createVector l
@@ -192,23 +196,23 @@
         joiner as tot ptr
     return r
   where joiner [] _ _ = return ()
-        joiner (V {dim = n, fptr = b} : cs) _ p = do
+        joiner (V {idim = n, fptr = b} : cs) _ p = do
             withForeignPtr b $ \pb -> copyArray p pb n
             joiner cs 0 (advancePtr p n)
 
 
 -- | transforms a complex vector into a real vector with alternating real and imaginary parts 
 asReal :: Vector (Complex Double) -> Vector Double
-asReal v = V { dim = 2*dim v, fptr =  castForeignPtr (fptr v) }
+asReal v = V { idim = 2*dim v, fptr =  castForeignPtr (fptr v) }
 
 -- | transforms a real vector into a complex vector with alternating real and imaginary parts
 asComplex :: Vector Double -> Vector (Complex Double)
-asComplex v = V { dim = dim v `div` 2, fptr =  castForeignPtr (fptr v) }
+asComplex v = V { idim = dim v `div` 2, fptr =  castForeignPtr (fptr v) }
 
 ----------------------------------------------------------------
 
 cloneVector :: Storable t => Vector t -> IO (Vector t)
-cloneVector (v@V {dim=n}) = do
+cloneVector (v@V {idim=n}) = do
         r <- createVector n
         let f _ s _ d =  copyArray d s n >> return 0
         app2 f vec v vec r "cloneVector"
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -27,8 +27,8 @@
     subMatrix, takeRows, dropRows, takeColumns, dropColumns,
     extractRows,
     ident, diag, diagRect, takeDiag,
-    liftMatrix, liftMatrix2,
-    format, dispf, disps,
+    liftMatrix, liftMatrix2, liftMatrix2Auto,
+    dispf, disps, dispcf, latexFormat, format,
     loadMatrix, saveMatrix, fromFile, fileDimensions,
     readMatrix, fromArray2D
 ) where
@@ -36,10 +36,11 @@
 import Data.Packed.Internal
 import qualified Data.Packed.ST as ST
 import Data.Packed.Vector
-import Data.List(transpose,intersperse)
 import Data.Array
 import System.Process(readProcess)
 import Text.Printf(printf)
+import Data.List(transpose,intersperse)
+import Data.Complex
 
 -- | creates a matrix from a vertical list of matrices
 joinVert :: Element t => [Matrix t] -> Matrix t
@@ -85,7 +86,7 @@
     rs = map (compatdim . map rows) ms
     cs = map (compatdim . map cols) (transpose ms)
     szs = sequence [rs,cs]
-    ms' = partit bc $ zipWith g szs (concat ms)
+    ms' = splitEvery bc $ zipWith g szs (concat ms)
 
     g [Just nr,Just nc] m
                 | nr == r && nc == c = m
@@ -232,30 +233,10 @@
           r = r1-r0+1
           c = c1-c0+1
 
-------------------------------------------------------
-{-
--- shows a Double with n digits after the decimal point    
-shf :: (RealFloat a) => Int -> a -> String     
-shf dec n | abs n < 1e-10 = "0."
-          | abs (n - (fromIntegral.round $ n)) < 1e-10 = show (round n) ++"."
-          | otherwise = showGFloat (Just dec) n ""    
--- shows a Complex Double as a pair, with n digits after the decimal point    
-shfc n z@ (a:+b) 
-    | magnitude z <1e-10 = "0."
-    | abs b < 1e-10 = shf n a
-    | abs a < 1e-10 = shf n b ++"i"
-    | b > 0         = shf n a ++"+"++shf n b ++"i"
-    | otherwise     = shf n a ++shf n b ++"i"         
 
--}
+-------------------------------------------------------------------
+-- display utilities
 
-dsp' :: String -> [[String]] -> String
-dsp' sep as = unlines . map unwords' $ transpose mtp where 
-    mt = transpose as
-    longs = map (maximum . map length) mt
-    mtp = zipWith (\a b -> map (pad a) b) longs mt
-    pad n str = replicate (n - length str) ' ' ++ str
-    unwords' = concat . intersperse sep
 
 {- | Creates a string from a matrix given a separator and a function to show each entry. Using
 this function the user can easily define any desired display function:
@@ -266,20 +247,7 @@
 
 -}
 format :: (Element t) => String -> (t -> String) -> Matrix t -> String
-format sep f m = dsp' sep . map (map f) . toLists $ m
-
-{-
-disp m f = putStrLn $ "matrix ("++show (rows m) ++"x"++ show (cols m) ++")\n"++format " | " f m
-
-dispR :: Int -> Matrix Double -> IO ()
-dispR d m = disp m (shf d)
-
-dispC :: Int -> Matrix (Complex Double) -> IO ()
-dispC d m = disp m (shfc d)
--}
-
--------------------------------------------------------------------
--- display utilities
+format sep f m = table sep . map (map f) . toLists $ m
 
 {- | Show a matrix with \"autoscaling\" and a given number of decimal places.
 
@@ -314,9 +282,7 @@
 
 formatFixed d x = format "  " (printf ("%."++show d++"f")) $ x
 
-isInt = all lookslikeInt . toList . flatten where
-    lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx
-        where shx = show x
+isInt = all lookslikeInt . toList . flatten
 
 formatScaled dec t = "E"++show o++"\n" ++ ss
     where ss = format " " (printf fmt. g) t
@@ -340,7 +306,43 @@
     . f . trans . reshape 1
     $ v
 
+-- | Tool to display matrices with latex syntax.
+latexFormat :: String -- ^ type of braces: \"matrix\", \"bmatrix\", \"pmatrix\", etc.
+            -> String -- ^ Formatted matrix, with elements separated by spaces and newlines
+            -> String
+latexFormat del tab = "\\begin{"++del++"}\n" ++ f tab ++ "\\end{"++del++"}"
+    where f = unlines . intersperse "\\\\" . map unwords . map (intersperse " & " . words) . tail . lines
 
+-- | Pretty print a complex number with at most n decimal digits.
+showComplex :: Int -> Complex Double -> String
+showComplex d (a:+b)
+    | isZero a && isZero b = "0"
+    | isZero b = sa
+    | isZero a && isOne b = s2++"i"
+    | isZero a = sb++"i"
+    | isOne b = sa++s3++"i"
+    | otherwise = sa++s1++sb++"i"
+  where
+    sa = shcr d a
+    sb = shcr d b
+    s1 = if b<0 then "" else "+"
+    s2 = if b<0 then "-" else ""
+    s3 = if b<0 then "-" else "+"
+
+shcr d a | lookslikeInt a = printf "%.0f" a
+         | otherwise      = printf ("%."++show d++"f") a
+
+
+lookslikeInt x = show (round x :: Int) ++".0" == shx || "-0.0" == shx
+   where shx = show x
+
+isZero x = show x `elem` ["0.0","-0.0"]
+isOne  x = show x `elem` ["1.0","-1.0"]
+
+-- | Pretty print a complex matrix with at most n decimal digits.
+dispcf :: Int -> Matrix (Complex Double) -> String
+dispcf d m = sdims m ++ "\n" ++ format "  " (showComplex d) m
+
 --------------------------------------------------------------------
 
 -- | reads a matrix from a string containing a table of numbers.
@@ -385,4 +387,29 @@
 
 -}
 repmat :: (Element t) => Matrix t -> Int -> Int -> Matrix t
-repmat m r c = fromBlocks $ partit c $ replicate (r*c) m
+repmat m r c = fromBlocks $ splitEvery c $ replicate (r*c) m
+
+-- | A version of 'liftMatrix2' which automatically adapt matrices with a single row or column to match the dimensions of the other matrix.
+liftMatrix2Auto :: (Element t, Element a, Element b)
+             => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
+liftMatrix2Auto f m1 m2 | compat' m1 m2 = lM f m1 m2
+                     | rows m1 == rows m2 && cols m2 == 1 = lM f m1 (repCols (cols m1) m2)
+                     | rows m1 == rows m2 && cols m1 == 1 = lM f (repCols (cols m2) m1) m2
+                     | cols m1 == cols m2 && rows m2 == 1 = lM f m1 (repRows (rows m1) m2)
+                     | cols m1 == cols m2 && cols m1 == 1 = lM f (repRows (rows m2) m1) m2
+                     | rows m1 == 1 && cols m2 == 1 = lM f (repRows (rows m2) m1) (repCols (cols m1) m2)
+                     | cols m1 == 1 && rows m2 == 1 = lM f (repCols (cols m2) m1) (repRows (rows m1) m2)
+                     | otherwise    = error $ "nonconformable matrices in liftMatrix2Auto: " ++ show (size m1) ++ ", " ++ show (size m2)
+
+size m = (rows m, cols m)
+
+lM f m1 m2 = reshape (max (cols m1) (cols m2)) (f (flatten m1) (flatten m2))
+
+repRows n x = fromRows (replicate n (flatten x))
+repCols n x = fromColumns (replicate n (flatten x))
+
+compat' :: Matrix a -> Matrix b -> Bool
+compat' m1 m2 = rows m1 == 1 && cols m1 == 1
+             || rows m2 == 1 && cols m2 == 1
+             || rows m1 == rows m2 && cols m1 == cols m2
+
diff --git a/lib/Numeric/GSL.hs b/lib/Numeric/GSL.hs
--- a/lib/Numeric/GSL.hs
+++ b/lib/Numeric/GSL.hs
@@ -19,6 +19,7 @@
 , module Numeric.GSL.Polynomials
 , module Numeric.GSL.Minimization
 , module Numeric.GSL.Root
+, module Numeric.GSL.ODE
 , module Numeric.GSL.Special
 , module Complex
 , setErrorHandlerOff
@@ -31,6 +32,7 @@
 import Numeric.GSL.Polynomials
 import Numeric.GSL.Minimization
 import Numeric.GSL.Root
+import Numeric.GSL.ODE
 import Complex
 
 
diff --git a/lib/Numeric/GSL/Internal.hs b/lib/Numeric/GSL/Internal.hs
--- a/lib/Numeric/GSL/Internal.hs
+++ b/lib/Numeric/GSL/Internal.hs
@@ -30,6 +30,9 @@
 foreign import ccall "wrapper"
     mkVecVecfun :: TVV -> IO (FunPtr TVV)
 
+foreign import ccall "wrapper"
+    mkDoubleVecVecfun :: (Double -> TVV) -> IO (FunPtr (Double -> TVV))
+
 aux_vTov :: (Vector Double -> Vector Double) -> TVV
 aux_vTov f n p nr r = g where
     V {fptr = pr} = f x
@@ -42,6 +45,9 @@
 
 foreign import ccall "wrapper"
     mkVecMatfun :: TVM -> IO (FunPtr TVM)
+
+foreign import ccall "wrapper"
+    mkDoubleVecMatfun :: (Double -> TVM) -> IO (FunPtr (Double -> TVM))
 
 aux_vTom :: (Vector Double -> Matrix Double) -> TVM
 aux_vTom f n p rr cr r = g where
diff --git a/lib/Numeric/GSL/ODE.hs b/lib/Numeric/GSL/ODE.hs
new file mode 100644
--- /dev/null
+++ b/lib/Numeric/GSL/ODE.hs
@@ -0,0 +1,111 @@
+{- |
+Module      :  Numeric.GSL.ODE
+Copyright   :  (c) Alberto Ruiz 2010
+License     :  GPL
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  uses ffi
+
+Solution of ordinary differential equation (ODE) initial value problems.
+
+<http://www.gnu.org/software/gsl/manual/html_node/Ordinary-Differential-Equations.html>
+
+A simple example:
+
+@import Numeric.GSL
+import Numeric.LinearAlgebra
+import Graphics.Plot
+
+xdot t [x,v] = [v, -0.95*x - 0.1*v]
+
+ts = linspace 100 (0,20)
+
+sol = odeSolve xdot [10,0] ts
+
+main = mplot (ts : toColumns sol)@
+
+-}
+-----------------------------------------------------------------------------
+
+module Numeric.GSL.ODE (
+    odeSolve, odeSolveV, ODEMethod(..)
+) where
+
+import Data.Packed.Internal
+import Data.Packed.Matrix
+import Foreign
+import Foreign.C.Types(CInt)
+import Numeric.GSL.Internal
+
+-------------------------------------------------------------------------
+
+-- | Stepping functions
+data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method.
+               | RK4 -- ^ 4th order (classical) Runge-Kutta. The error estimate is obtained by halving the step-size. For more efficient estimate of the error, use 'RKf45'.
+               | RKf45 -- ^ Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator.
+               | RKck -- ^ Embedded Runge-Kutta Cash-Karp (4, 5) method.
+               | RK8pd -- ^ Embedded Runge-Kutta Prince-Dormand (8,9) method.
+               | RK2imp -- ^ Implicit 2nd order Runge-Kutta at Gaussian points.
+               | RK4imp -- ^ Implicit 4th order Runge-Kutta at Gaussian points.
+               | BSimp -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. This algorithm requires the Jacobian.
+               | Gear1 -- ^ M=1 implicit Gear method.
+               | Gear2 -- ^ M=2 implicit Gear method.
+               deriving (Enum,Eq,Show,Bounded)
+
+-- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists.
+odeSolve
+    :: (Double -> [Double] -> [Double])        -- ^ xdot(t,x)
+    -> [Double]        -- ^ initial conditions
+    -> Vector Double   -- ^ desired solution times
+    -> Matrix Double   -- ^ solution
+odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) Nothing (fromList xi) ts
+    where hi = (ts@>1 - ts@>0)/100
+          epsAbs = 1.49012e-08
+          epsRel = 1.49012e-08
+          l2v f = \t -> fromList  . f t . toList
+          l2m f = \t -> fromLists . f t . toList
+
+-- | Evolution of the system with adaptive step-size control.
+odeSolveV
+    :: ODEMethod
+    -> Double -- ^ initial step size
+    -> Double -- ^ absolute tolerance for the state vector
+    -> Double -- ^ relative tolerance for the state vector
+    -> (Double -> Vector Double -> Vector Double)   -- ^ xdot(t,x)
+    -> Maybe (Double -> Vector Double -> Matrix Double)   -- ^ optional jacobian
+    -> Vector Double     -- ^ initial conditions
+    -> Vector Double     -- ^ desired solution times
+    -> Matrix Double     -- ^ solution
+odeSolveV method h epsAbs epsRel f mbjac xiv ts = unsafePerformIO $ do
+    let n   = dim xiv
+    fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t))
+    jp <- case mbjac of
+        Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t))
+        Nothing  -> return nullFunPtr
+    sol <- withVector xiv $ \xiv' ->
+            withVector (checkTimes ts) $ \ts' ->
+             createMIO (dim ts) n
+              (ode_c (fi (fromEnum method)) h epsAbs epsRel fp jp // xiv' // ts' )
+              "ode"
+    freeHaskellFunPtr fp
+    return sol
+
+foreign import ccall "ode"
+    ode_c :: CInt -> Double -> Double -> Double -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVM
+
+-------------------------------------------------------
+
+checkdim1 n v
+    | dim v == n = v
+    | otherwise = error $ "Error: "++ show n
+                        ++ " components expected in the result of the function supplied to odeSolve"
+
+checkdim2 n m
+    | rows m == n && cols m == n = m
+    | otherwise = error $ "Error: "++ show n ++ "x" ++ show n
+                        ++ " Jacobian expected in odeSolve"
+
+checkTimes ts | dim ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts
+              | otherwise = error "odeSolve requires increasing times"
+    where ts' = toList ts
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -22,6 +22,7 @@
 #include <gsl/gsl_complex_math.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
+#include <gsl/gsl_odeiv.h>
 #include <string.h>
 #include <stdio.h>
 
@@ -507,7 +508,7 @@
 
 
 void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) {
-    Tfdf * fdf = ((Tfdf*) pars);  
+    Tfdf * fdf = ((Tfdf*) pars);
     double* p = (double*)calloc(x->size,sizeof(double));
     double* q = (double*)calloc(g->size,sizeof(double));
     int k;
@@ -797,3 +798,89 @@
     }
 }
 #undef RAN
+
+//////////////////////////////////////////////////////
+//                        ODE
+
+typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode;
+
+int odefunc (double t, const double y[], double f[], void *params) { 
+    Tode * P = (Tode*) params;
+    (P->f)(t,P->n,y,P->n,f);
+    return GSL_SUCCESS;
+}
+
+int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) {
+     Tode * P = ((Tode*) params);
+     (P->j)(t,P->n,y,P->n,P->n,dfdy);
+     int j;
+     for (j=0; j< P->n; j++)
+        dfdt[j] = 0.0;
+     return GSL_SUCCESS;
+}
+
+
+int ode(int method, double h, double eps_abs, double eps_rel,
+        int f(double, int, const double*, int, double*),
+        int jac(double, int, const double*, int, int, double*),
+        KRVEC(xi), KRVEC(ts), RMAT(sol)) {
+
+    const gsl_odeiv_step_type * T;
+
+    switch(method) {
+        case 0 : {T = gsl_odeiv_step_rk2; break; }
+        case 1 : {T = gsl_odeiv_step_rk4; break; }
+        case 2 : {T = gsl_odeiv_step_rkf45; break; }
+        case 3 : {T = gsl_odeiv_step_rkck; break; }
+        case 4 : {T = gsl_odeiv_step_rk8pd; break; }
+        case 5 : {T = gsl_odeiv_step_rk2imp; break; }
+        case 6 : {T = gsl_odeiv_step_rk4imp; break; }
+        case 7 : {T = gsl_odeiv_step_bsimp; break; }
+        case 8 : {T = gsl_odeiv_step_gear1; break; }
+        case 9 : {T = gsl_odeiv_step_gear2; break; }
+        default: ERROR(BAD_CODE);
+    }
+
+
+    gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin);
+    gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel);
+    gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin);
+
+    Tode P;
+    P.f = f;
+    P.j = jac;
+    P.n = xin;
+
+    gsl_odeiv_system sys = {odefunc, odejac, xin, &P};
+
+    double t = tsp[0];
+
+    double* y = (double*)calloc(xin,sizeof(double));
+    int i,j;
+    for(i=0; i< xin; i++) {
+        y[i] = xip[i];
+        solp[i] = xip[i];
+    }
+
+       for (i = 1; i < tsn ; i++)
+         {
+           double ti = tsp[i];
+           while (t < ti)
+             {
+               gsl_odeiv_evolve_apply (e, c, s,
+                                       &sys,
+                                       &t, ti, &h,
+                                       y);
+               // if (h < hmin) h = hmin;
+             }
+           for(j=0; j<xin; j++) {
+               solp[i*xin + j] = y[j];
+           }
+         }
+
+    free(y);
+    gsl_odeiv_evolve_free (e);
+    gsl_odeiv_control_free (c);
+    gsl_odeiv_step_free (s);
+    return 0;
+}
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -729,7 +729,7 @@
 -}
 kronecker :: (Field t) => Matrix t -> Matrix t -> Matrix t
 kronecker a b = fromBlocks
-              . partit (cols a)
+              . splitEvery (cols a)
               . map (reshape (cols b))
               . toRows
               $ flatten a `outer` flatten b
diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs
--- a/lib/Numeric/LinearAlgebra/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Instances.hs
@@ -71,15 +71,7 @@
     | dim y == 1 = f3 x (y@>0)
     | otherwise = f2 x y
 
-liftMatrix2' :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
-liftMatrix2' f m1 m2 | compat' m1 m2 = reshape (max (cols m1) (cols m2)) (f (flatten m1) (flatten m2))
-                     | otherwise    = error "nonconformant matrices in liftMatrix2'"
 
-compat' :: Matrix a -> Matrix b -> Bool
-compat' m1 m2 = rows m1 == 1 && cols m1 == 1
-             || rows m2 == 1 && cols m2 == 1
-             || rows m1 == rows m2 && cols m1 == cols m2
-
 instance Linear Vector a => Eq (Vector a) where
     (==) = equal
 
@@ -103,10 +95,10 @@
     (==) = equal
 
 instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where
-    (+) = liftMatrix2' (+)
-    (-) = liftMatrix2' (-)
+    (+) = liftMatrix2Auto (+)
+    (-) = liftMatrix2Auto (-)
     negate = liftMatrix negate
-    (*) = liftMatrix2' (*)
+    (*) = liftMatrix2Auto (*)
     signum = liftMatrix signum
     abs = liftMatrix abs
     fromInteger = (1><1) . return . fromInteger
@@ -123,7 +115,7 @@
 
 instance (Linear Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where
     fromRational n = (1><1) [fromRational n]
-    (/) = liftMatrix2' (/)
+    (/) = liftMatrix2Auto (/)
 
 ---------------------------------------------------------
 
@@ -184,14 +176,14 @@
     atanh = liftMatrix atanh
     exp   = liftMatrix exp
     log   = liftMatrix log
-    (**)  = liftMatrix2' (**)
+    (**)  = liftMatrix2Auto (**)
     sqrt  = liftMatrix sqrt
     pi    = (1><1) [pi]
 
 ---------------------------------------------------------------
 
 instance (Storable a, Num (Vector a)) => Monoid (Vector a) where
-    mempty = 0 { dim = 0 }
+    mempty = 0 { idim = 0 }
     mappend a b = mconcat [a,b]
     mconcat = j . filter ((>0).dim)
         where j [] = mempty
@@ -204,3 +196,4 @@
 --
 -- instance (NFData a, Element a) => NFData (Matrix a) where
 --     rnf = rnf . flatten
+
diff --git a/lib/Numeric/LinearAlgebra/Interface.hs b/lib/Numeric/LinearAlgebra/Interface.hs
--- a/lib/Numeric/LinearAlgebra/Interface.hs
+++ b/lib/Numeric/LinearAlgebra/Interface.hs
@@ -72,6 +72,9 @@
 
 ------------------------------------------------
 
+{-# DEPRECATED (<|>) "define operator a & b = fromBlocks[[a,b]] and use asRow/asColumn to join vectors" #-}
+{-# DEPRECATED (<->) "define operator a // b = fromBlocks[[a],[b]] and use asRow/asColumn to join vectors" #-}
+
 class Joinable a b where
     joinH :: Element t => a t -> b t -> Matrix t
     joinV :: Element t => a t -> b t -> Matrix t
@@ -108,3 +111,4 @@
 -- | Vertical concatenation of matrices and vectors.
 (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
 a <-> b = joinV a b
+
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -135,6 +135,14 @@
 
 ---------------------------------------------------------------------
 
+odeTest = utest "ode" (last (toLists sol) ~~ [-1.7588880332411019, 8.364348908711941e-2])
+    where sol = odeSolveV RK8pd 1E-6 1E-6 0 (l2v $ vanderpol 10) Nothing (fromList [1,0]) ts
+          ts = linspace 101 (0,100)
+          l2v f = \t -> fromList  . f t . toList
+          vanderpol mu _t [x,y] = [y, -x + mu * y * (1-x^2) ]
+
+---------------------------------------------------------------------
+
 randomTestGaussian = c :~1~: snd (meanCov dat) where
     a = (3><3) [1,2,3,
                 2,4,0,
@@ -280,6 +288,7 @@
         , utest "rank" $  rank ((2><3)[1,0,0,1,6*eps,0]) == 1
                        && rank ((2><3)[1,0,0,1,7*eps,0]) == 2
         , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM)
+        , odeTest
         ]
     return ()
 
@@ -366,3 +375,4 @@
     time "full svd         3000x500" (fv $ svd a)
     time "singular values 1000x1000" (singularValues b)
     time "full svd        1000x1000" (fv $ svd b)
+
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
--- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -148,7 +148,7 @@
            && orthonormal u && orthonormal v
            && (dim s == r || r == 0 && dim s == 1)
     where (u,s,v) = compactSVD m
-          m = m' <-> m'
+          m = fromBlocks [[m'],[m']]
           r = rank m'
 
 svdProp5a m = and (map (s1|~|) [s2,s3,s4,s5,s6]) where