diff --git a/README-quick-start.md b/README-quick-start.md
new file mode 100644
index 000000000..2adc04d70
--- /dev/null
+++ b/README-quick-start.md
@@ -0,0 +1,606 @@
+# Quickstart tutorial
+
+If you are familiar with Python Numpy, do check out this [For Numpy User Doc](https://docs.rs/ndarray/0.13.0/ndarray/doc/ndarray_for_numpy_users/index.html)
+after you go through this tutorial. 
+
+You can use [play.integer32.com](https://play.integer32.com/) to immediately try out the examples.
+
+## The Basics
+
+Just create your first 2x3 floating-point ndarray 
+```rust
+use ndarray::prelude::*;
+
+fn main() {
+    let a = array![
+                [1.,2.,3.], 
+                [4.,5.,6.],
+            ]; 
+    assert_eq!(a.ndim(), 2);         // get the number of dimensions of array a
+    assert_eq!(a.len(), 6);          // get the number of elements in array a
+    assert_eq!(a.shape(), [2, 3]);   // get the shape of array a
+    assert_eq!(a.is_empty(), false); // check if the array has zero elements
+
+    println!("{:?}", a);
+}
+```
+This code will create a simple array and output to stdout:
+```
+[[1.0, 2.0, 3.0],
+ [4.0, 5.0, 6.0]], shape=[2, 3], strides=[3, 1], layout=C (0x1), const ndim=2
+```
+
+## Array Creation
+
+### Element type and dimensionality
+
+Now let's create more arrays. How about try make a zero array with dimension of (3, 2, 4)?
+
+```rust
+use ndarray::prelude::*;
+use ndarray::Array;
+fn main() {
+    let a = Array::zeros((3, 2, 4).f());
+    println!("{:?}", a);
+}
+```
+gives
+```
+|    let a = Array::zeros((3, 2, 4).f());
+|        -   ^^^^^^^^^^^^ cannot infer type for type parameter `A`
+```
+Note that the compiler needs to infer the element type and dimensionality from context. In this 
+case the compiler failed to do that. Now we give it the type and let it infer dimensionality
+
+```rust
+use ndarray::prelude::*;
+use ndarray::Array;
+fn main() {
+  let a = Array::<f64, _>::zeros((3, 2, 4).f());
+  println!("{:?}", a);
+}
+```
+and now it works:
+```
+[[[0.0, 0.0, 0.0, 0.0],
+  [0.0, 0.0, 0.0, 0.0]],
+
+ [[0.0, 0.0, 0.0, 0.0],
+  [0.0, 0.0, 0.0, 0.0]],
+
+ [[0.0, 0.0, 0.0, 0.0],
+  [0.0, 0.0, 0.0, 0.0]]], shape=[3, 2, 4], strides=[1, 3, 6], layout=F (0x2), const ndim=3
+```
+
+We can also specify its dimensionality
+
+```rust
+use ndarray::prelude::*;
+use ndarray::{Array, Ix3};
+fn main() {
+  let a = Array::<f64, Ix3>::zeros((3, 2, 4).f());
+  println!("{:?}", a);
+}
+```
+`Ix3` stands for 3D array.
+
+And now we are type checked. Try change the code above to `Array::<f64, Ix3>::zeros((3, 2, 4, 5).f());`
+and compile, see what happens.
+
+### How about create array of different type and having different initial values?
+
+The [`from_elm`](http://docs.rs/ndarray/latest/ndarray/struct.ArrayBase.html#method.from_elem) method can be handy here:
+
+```rust
+use ndarray::{Array, Ix3};
+fn main() {
+  let a = Array::<bool, Ix3>::from_elem((3, 2, 4), false);
+  println!("{:?}", a);
+}
+```
+
+### Some common create helper functions
+`linspace` - Create a 1-D array with 21 elements with values 0., …, 5.
+```rust
+use ndarray::prelude::*;
+use ndarray::{Array, Ix3};
+fn main() {
+  let a = Array::<f64, _>::linspace(0., 5., 11);
+  println!("{:?}", a);
+}
+```
+The output is:
+```
+[0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0], shape=[11], strides=[1], layout=C | F (0x3), const ndim=1
+```
+
+And there are also `range`, `logspace`, `ones`, `eye` and so on you can choose to use.
+
+## Basic operations
+
+```rust
+use ndarray::prelude::*;
+use ndarray::Array;
+use std::f64::INFINITY as inf;
+
+fn main() {
+    let a = array![
+                [10.,20.,30., 40.,], 
+            ];
+    let b = Array::range(0., 4., 1.);  // [0., 1., 2., 3, ]
+
+    assert_eq!(&a + &b, array![[10., 21., 32., 43.,]]);  // Allocates a new array. Note the explicit `&`.
+    assert_eq!(&a - &b, array![[10., 19., 28., 37.,]]);
+    assert_eq!(&a * &b, array![[0., 20., 60., 120.,]]);
+    assert_eq!(&a / &b, array![[inf, 20., 15., 13.333333333333334,]]);
+}
+```
+
+Try remove all the `&` sign in front of `a` and `b`, does it still compile? Why?
+
+Note that
+* `&A @ &A` produces a new `Array`
+* `B @ A` consumes `B`, updates it with the result, and returns it
+* `B @ &A` consumes `B`, updates it with the result, and returns it
+* `C @= &A` performs an arithmetic operation in place
+
+For more info checkout https://docs.rs/ndarray/latest/ndarray/struct.ArrayBase.html#arithmetic-operations
+
+Some operations have `_axis` appended to the function name: they generally take in a parameter of type `Axis` as one of their inputs,
+such as `sum_axis`:
+
+```rust
+use ndarray::{aview0, aview1, arr2, Axis};
+
+fn main() {
+    let a = arr2(&[[1., 2., 3.],
+                   [4., 5., 6.]]);
+    assert!(
+        a.sum_axis(Axis(0)) == aview1(&[5., 7., 9.]) &&
+        a.sum_axis(Axis(1)) == aview1(&[6., 15.]) &&
+
+        a.sum_axis(Axis(0)).sum_axis(Axis(0)) == aview0(&21.) &&
+        a.sum_axis(Axis(0)).sum_axis(Axis(0)) == aview0(&a.sum())
+    );
+}
+```
+
+### Matrix product
+
+```rust
+use ndarray::prelude::*;
+use ndarray::Array;
+
+fn main() {
+    let a = array![
+                [10.,20.,30., 40.,], 
+            ];
+    let b = Array::range(0., 4., 1.);     // b = [0., 1., 2., 3, ]
+    println!("a shape {:?}", &a.shape());
+    println!("b shape {:?}", &b.shape());
+    
+    let b = b.into_shape((4,1)).unwrap(); // reshape b to shape [4, 1]
+    println!("b shape {:?}", &b.shape());
+    
+    println!("{}", a.dot(&b));            // [1, 4] x [4, 1] -> [1, 1] 
+    println!("{}", a.t().dot(&b.t()));    // [4, 1] x [1, 4] -> [4, 4]
+}
+```
+The output is:
+```
+a shape [1, 4]
+b shape [4]
+b shape after reshape [4, 1]
+[[200]]
+[[0, 10, 20, 30],
+ [0, 20, 40, 60],
+ [0, 30, 60, 90],
+ [0, 40, 80, 120]]
+```
+
+## Indexing, Slicing and Iterating
+One-dimensional arrays can be indexed, sliced and iterated over, much like `numpy` arrays
+
+```rust
+use ndarray::prelude::*;
+use ndarray::Array;
+
+fn main() {
+    let a = Array::range(0., 10., 1.);
+
+    let mut a = a.mapv(|a: f64| a.powi(3));  // numpy equivlant of `a ** 3`; https://doc.rust-lang.org/nightly/std/primitive.f64.html#method.powi
+
+    println!("{}", a);
+
+    println!("{}", a[[2]]);
+    println!("{}", a.slice(s![2]));
+
+    println!("{}", a.slice(s![2..5]));
+
+    a.slice_mut(s![..6;2]).fill(1000.);  // numpy equivlant of `a[:6:2] = 1000`
+    println!("{}", a);
+
+    for i in a.iter() {
+        print!("{}, ", i.powf(1./3.))
+    }
+}
+```
+The output is:
+```
+[0, 1, 8, 27, 64, 125, 216, 343, 512, 729]
+8
+8
+[8, 27, 64]
+[1000, 1, 1000, 27, 1000, 125, 216, 343, 512, 729]
+9.999999999999998, 1, 9.999999999999998, 3, 9.999999999999998, 4.999999999999999, 5.999999999999999, 6.999999999999999, 7.999999999999999, 8.999999999999998,
+```
+
+For more info about iteration see [Loops, Producers, and Iterators](https://docs.rs/ndarray/0.13.0/ndarray/struct.ArrayBase.html#loops-producers-and-iterators)
+
+Let's try a 3D array with elements of type `isize`. This is how you index it: 
+```rust
+use ndarray::prelude::*;
+
+fn main() {
+    let a = array![
+                    [[  0,  1,  2],         // a 3D array  2 x 2 x 3
+                     [ 10, 12, 13]],
+
+                    [[100,101,102],
+                     [110,112,113]]
+                ];
+
+    let a = a.mapv(|a: isize| a.pow(1));  // numpy equivlant of `a ** 1`; 
+                                          // This line does nothing but illustrate mapv with isize type 
+    println!("a -> \n{}\n", a);
+
+    println!("`a.slice(s![1, .., ..])` -> \n{}\n", a.slice(s![1, .., ..]));
+
+    println!("`a.slice(s![.., .., 2])` -> \n{}\n", a.slice(s![.., .., 2]));
+
+    println!("`a.slice(s![.., 1, 0..2])` -> \n{}\n", a.slice(s![.., 1, 0..2]));
+
+    println!("`a.iter()` ->");
+    for i in a.iter() {
+        print!("{}, ", i)  // flat out to every element
+    }
+
+    println!("\n\n`a.outer_iter()` ->");
+    for i in a.outer_iter() {
+        print!("row: {}, \n", i)  // iterate through first dimension
+    }
+}
+```
+The output is:
+```
+a -> 
+[[[0, 1, 2],
+  [10, 12, 13]],
+
+ [[100, 101, 102],
+  [110, 112, 113]]]
+
+`a.slice(s![1, .., ..])` -> 
+[[100, 101, 102],
+ [110, 112, 113]]
+
+`a.slice(s![.., .., 2])` -> 
+[[2, 13],
+ [102, 113]]
+
+`a.slice(s![.., 1, 0..2])` -> 
+[[10, 12],
+ [110, 112]]
+
+`a.iter()` ->
+0, 1, 2, 10, 12, 13, 100, 101, 102, 110, 112, 113, 
+
+`a.outer_iter()` ->
+row: [[0, 1, 2],
+ [10, 12, 13]], 
+row: [[100, 101, 102],
+ [110, 112, 113]], 
+```
+
+## Shape Manipulation
+
+### Changing the shape of an array
+The shape of an array can be changed with `into_shape` method.
+
+````rust
+use ndarray::prelude::*;
+use ndarray::Array;
+use std::iter::FromIterator;
+// use ndarray_rand::RandomExt;
+// use ndarray_rand::rand_distr::Uniform;
+
+fn main() {
+    // Or you may use ndarray_rand crate to generate random arrays
+    // let a = Array::random((2, 5), Uniform::new(0., 10.));
+    
+    let a = array![
+        [3., 7., 3., 4.],
+        [1., 4., 2., 2.],
+        [7., 2., 4., 9.]];
+        
+    println!("a = \n{:?}\n", a);
+    
+    // use trait FromIterator to flatten a matrix to a vector
+    let b = Array::from_iter(a.iter());
+    println!("b = \n{:?}\n", b);
+    
+    let c = b.into_shape([6, 2]).unwrap(); // consume b and generate c with new shape
+    println!("c = \n{:?}", c);
+}
+````
+The output is:
+```
+a = 
+[[3.0, 7.0, 3.0, 4.0],
+ [1.0, 4.0, 2.0, 2.0],
+ [7.0, 2.0, 4.0, 9.0]], shape=[3, 4], strides=[4, 1], layout=C (0x1), const ndim=2
+
+b = 
+[3.0, 7.0, 3.0, 4.0, 1.0, 4.0, 2.0, 2.0, 7.0, 2.0, 4.0, 9.0], shape=[12], strides=[1], layout=C | F (0x3), const ndim=1
+
+c = 
+[[3.0, 7.0],
+ [3.0, 4.0],
+ [1.0, 4.0],
+ [2.0, 2.0],
+ [7.0, 2.0],
+ [4.0, 9.0]], shape=[6, 2], strides=[2, 1], layout=C (0x1), const ndim=2
+```
+
+### Stacking together different arrays
+
+Macro `stack!` is helpful for stacking arrays: 
+```rust
+use ndarray::prelude::*;
+use ndarray::{Array, Axis, stack};
+
+fn main() {
+
+    let a = array![
+        [9., 7.],
+        [5., 2.]];
+    
+    let b = array![
+        [1., 9.],
+        [5., 1.]];
+        
+    println!("a vstack b = \n{:?}\n", stack![Axis(0), a, b]);
+    
+    println!("a hstack b = \n{:?}\n", stack![Axis(1), a, b]);
+}
+```
+The output is:
+```
+a vstack b = 
+[[9.0, 7.0],
+ [5.0, 2.0],
+ [1.0, 9.0],
+ [5.0, 1.0]], shape=[4, 2], strides=[2, 1], layout=C (0x1), const ndim=2
+
+a hstack b = 
+[[9.0, 7.0, 1.0, 9.0],
+ [5.0, 2.0, 5.0, 1.0]], shape=[2, 4], strides=[4, 1], layout=C (0x1), const ndim=2
+```
+
+### Splitting one array into several smaller ones
+
+More to see here [ArrayView::split_at](https://docs.rs/ndarray/latest/ndarray/type.ArrayView.html#method.split_at)
+```rust
+use ndarray::prelude::*;
+use ndarray::Axis;
+
+fn main() {
+
+    let a = array![
+        [6., 7., 6., 9., 0., 5., 4., 0., 6., 8., 5., 2.],
+        [8., 5., 5., 7., 1., 8., 6., 7., 1., 8., 1., 0.]];
+    
+    let (s1, s2) = a.view().split_at(Axis(0), 1);
+    println!("Split a from Axis(0), at index 1:");
+    println!("s1  = \n{}", s1);
+    println!("s2  = \n{}\n", s2);
+    
+    
+    let (s1, s2) = a.view().split_at(Axis(1), 4);
+    println!("Split a from Axis(1), at index 4:");
+    println!("s1  = \n{}", s1);
+    println!("s2  = \n{}\n", s2);
+}
+```
+The output is:
+```
+Split a from Axis(0), at index 1:
+s1  = 
+[[6, 7, 6, 9, 0, 5, 4, 0, 6, 8, 5, 2]]
+s2  = 
+[[8, 5, 5, 7, 1, 8, 6, 7, 1, 8, 1, 0]]
+
+Split a from Axis(1), at index 4:
+s1  = 
+[[6, 7, 6, 9],
+ [8, 5, 5, 7]]
+s2  = 
+[[0, 5, 4, 0, 6, 8, 5, 2],
+ [1, 8, 6, 7, 1, 8, 1, 0]]
+
+```
+
+## Copies and Views
+### View, Ref or Shallow Copy
+As in Rust we have owner ship, so we cannot simply 
+update an element of an array while we have a 
+shared view of it. This will help us write more
+robust code.
+
+```rust
+use ndarray::prelude::*;
+use ndarray::{Array, Axis};
+
+fn main() {
+
+    let mut a = Array::range(0., 12., 1.).into_shape([3 ,4]).unwrap();
+    println!("a = \n{}\n", a);
+    
+    {
+        let (s1, s2) = a.view().split_at(Axis(1), 2);
+        
+        // with s as a view sharing the ref of a, we cannot update a here
+        // a.slice_mut(s![1, 1]).fill(1234.);
+        
+        println!("Split a from Axis(0), at index 1:");
+        println!("s1  = \n{}", s1);
+        println!("s2  = \n{}\n", s2);
+    }
+    
+    // now we can update a again here, as views of s1, s2 are dropped already
+    a.slice_mut(s![1, 1]).fill(1234.);
+    
+    let (s1, s2) = a.view().split_at(Axis(1), 2);
+    println!("Split a from Axis(0), at index 1:");
+    println!("s1  = \n{}", s1);
+    println!("s2  = \n{}\n", s2);
+}
+```
+The output is:
+```
+a = 
+[[0, 1, 2, 3],
+ [4, 5, 6, 7],
+ [8, 9, 10, 11]]
+
+Split a from Axis(0), at index 1:
+s1  = 
+[[0, 1],
+ [4, 5],
+ [8, 9]]
+s2  = 
+[[2, 3],
+ [6, 7],
+ [10, 11]]
+
+Split a from Axis(0), at index 1:
+s1  = 
+[[0, 1],
+ [4, 1234],
+ [8, 9]]
+s2  = 
+[[2, 3],
+ [6, 7],
+ [10, 11]]
+```
+
+### Deep Copy
+As the usual way in Rust, a `clone()` call will
+make a copy of your array:
+```rust
+use ndarray::prelude::*;
+use ndarray::Array;
+
+fn main() {
+
+    let mut a = Array::range(0., 4., 1.).into_shape([2 ,2]).unwrap();
+    let b = a.clone();
+    
+    println!("a = \n{}\n", a);
+    println!("b clone of a = \n{}\n", a);
+    
+    a.slice_mut(s![1, 1]).fill(1234.);
+    
+    println!("a updated...");
+    println!("a = \n{}\n", a);
+    println!("b clone of a = \n{}\n", b);
+}
+```
+
+The output is:
+```
+a = 
+[[0, 1],
+ [2, 3]]
+
+b clone of a = 
+[[0, 1],
+ [2, 3]]
+
+a updated...
+a = 
+[[0, 1],
+ [2, 1234]]
+
+b clone of a = 
+[[0, 1],
+ [2, 3]]
+```
+
+Noticing that using `clone()` (or cloning) an `Array` type also copies the array's elements. It creates an independently owned array of the same type.
+
+Cloning an `ArrayView` does not clone or copy the underlying elements - it just clones the view reference (as it happens in Rust when cloning a `&` reference).
+
+## Broadcasting
+
+Arrays support limited broadcasting, where arithmetic operations with array operands of different sizes can be carried out by repeating the elements of the smaller dimension array. 
+
+```rust
+use ndarray::prelude::*;
+
+fn main() {
+    let a = array![
+        [1., 1.], 
+        [1., 2.], 
+        [0., 3.], 
+        [0., 4.]];
+
+    let b = array![[0., 1.]];
+
+    let c = array![
+        [1., 2.], 
+        [1., 3.], 
+        [0., 4.], 
+        [0., 5.]];
+    
+    // We can add because the shapes are compatible even if not equal.
+    // The `b` array is shape 1 × 2 but acts like a 4 × 2 array.
+    assert!(c == a + b);
+}
+```
+
+See [.broadcast()](https://docs.rs/ndarray/latest/ndarray/struct.ArrayBase.html#method.broadcast) for a more detailed description.
+
+And there is a short example of it:
+```rust
+use ndarray::prelude::*;
+
+fn main() {
+    let a = array![
+        [1., 2.],
+        [3., 4.],
+    ];
+    
+    let b =  a.broadcast((3, 2, 2)).unwrap();
+    println!("shape of a is {:?}", a.shape());
+    println!("a is broadcased to 3x2x2 = \n{}", b);
+}
+```
+The output is:
+```
+shape of a is [2, 2]
+a is broadcased to 3x2x2 = 
+[[[1, 2],
+  [3, 4]],
+
+ [[1, 2],
+  [3, 4]],
+
+ [[1, 2],
+  [3, 4]]]
+```
+
+## Want to learn more?
+Please checkout these docs for more information
+* [`ArrayBase` doc page](https://docs.rs/ndarray/latest/ndarray/struct.ArrayBase.html)
+* [`ndarray` for `numpy` user doc page](https://docs.rs/ndarray/latest/ndarray/doc/ndarray_for_numpy_users/index.html)
diff --git a/README.rst b/README.rst
index f1a3524ad..9e2580a77 100644
--- a/README.rst
+++ b/README.rst
@@ -5,6 +5,7 @@ The ``ndarray`` crate provides an *n*-dimensional container for general elements
 and for numerics.
 
 Please read the `API documentation on docs.rs`__
+or take a look at the `quickstart tutorial <./README-quick-start.md>`_.
 
 __ https://docs.rs/ndarray/