diff --git a/docs/examples/portfolio.rst b/docs/examples/portfolio.rst
new file mode 100644
index 0000000..60649b6
--- /dev/null
+++ b/docs/examples/portfolio.rst
@@ -0,0 +1,165 @@
+Portfolio optimization
+======================
+
+
+Portfolio optimization seeks to allocate assets in a way that maximizes the risk adjusted return,
+
+
+.. math::
+  \begin{array}{ll}
+    \mbox{maximize} & \mu^T x - \gamma \left( x^T \Sigma x \right) \\
+    \mbox{subject to} & \boldsymbol{1}^T x = 1 \\
+                      & x \ge 0
+  \end{array}
+
+
+where :math:`x \in \mathbf{R}^{n}` represents the portfolio, :math:`\mu \in \mathbf{R}^{n}` the vector of expected returns, :math:`\gamma > 0` the risk aversion parameter, and :math:`\Sigma \in \mathbf{S}^{n}_{+}` the risk model covariance matrix.
+The risk model is usually assumed to be the sum of a diagonal and a rank :math:`k < n` matrix,
+
+
+.. math::
+  \Sigma = F F^T + D,
+
+
+where :math:`F \in \mathbf{R}^{n \times k}` is the factor loading matrix and :math:`D \in \mathbf{S}^{n}_{+}` is a diagonal matrix describing the asset-specific risk.
+The resulting problem has the following equivalent form,
+
+.. math::
+  \begin{array}{ll}
+    \mbox{minimize} & \frac{1}{2} x^T D x + \frac{1}{2} y^T y - \frac{1}{2\gamma}\mu^T x \\
+    \mbox{subject to} & y = F^T x \\
+                      & \boldsymbol{1}^T x = 1 \\
+                      & x \ge 0
+  \end{array}
+
+
+
+Python
+------
+
+.. code:: python
+
+    import osqp
+    import numpy as np
+    import scipy as sp
+    from scipy import sparse
+
+    # Generate problem data
+    sp.random.seed(1)
+    n = 100
+    k = 10
+    F = sparse.random(n, k, density=0.7, format='csc')
+    D = sparse.diags(np.random.rand(n) * np.sqrt(k), format='csc')
+    mu = np.random.randn(n)
+    gamma = 1
+
+    # OSQP data
+    P = sparse.block_diag([D, sparse.eye(k)], format='csc')
+    q = np.hstack([-mu / (2*gamma), np.zeros(k)])
+    A = sparse.vstack([
+            sparse.hstack([F.T, -sparse.eye(k)]),
+            sparse.hstack([sparse.csc_matrix(np.ones((1, n))), sparse.csc_matrix((1, k))]),
+            sparse.hstack((sparse.eye(n), sparse.csc_matrix((n, k))))
+        ], format='csc')
+    l = np.hstack([np.zeros(k), 1., np.zeros(n)])
+    u = np.hstack([np.zeros(k), 1., np.ones(n)])
+
+    # Create an OSQP object
+    prob = osqp.OSQP()
+
+    # Setup workspace
+    prob.setup(P, q, A, l, u)
+
+    # Solve problem
+    res = prob.solve()
+
+
+
+Matlab
+------
+
+.. code:: matlab
+
+    % Generate problem data
+    rng(1)
+    n = 100;
+    k = 10;
+    F = sprandn(n, k, 0.7);
+    D = sparse(diag( sqrt(k)*rand(n,1) ));
+    mu = randn(n, 1);
+    gamma = 1;
+
+    % OSQP data
+    P = blkdiag(D, speye(k));
+    q = [-mu/(2*gamma); zeros(k, 1)];
+    A = [F', -speye(k);
+         ones(1, n), zeros(1, k);
+         speye(n), sparse(n, k)];
+    l = [zeros(k, 1); 1; zeros(n, 1)];
+    u = [zeros(k, 1); 1; ones(n, 1)];
+
+    % Create an OSQP object
+    prob = osqp;
+
+    % Setup workspace
+    prob.setup(P, q, A, l, u);
+
+    % Solve problem
+    res = prob.solve();
+
+
+
+CVXPY
+-----
+
+.. code:: python
+
+    from cvxpy import *
+    import numpy as np
+    import scipy as sp
+    from scipy import sparse
+
+    # Generate problem data
+    sp.random.seed(1)
+    n = 100
+    k = 10
+    F = sparse.random(n, k, density=0.7, format='csc')
+    D = sparse.diags(np.random.rand(n) * np.sqrt(k), format='csc')
+    mu = np.random.randn(n)
+    gamma = 1
+    Sigma = F*F.T + D
+
+    # Define problem
+    x = Variable(n)
+    objective = mu.T*x - gamma*quad_form(x, Sigma)
+    constraints = [sum(x) == 1, x >= 0]
+
+    # Solve with OSQP
+    Problem(Maximize(objective), constraints).solve(solver=OSQP)
+
+
+
+YALMIP
+------
+
+.. code:: matlab
+
+    % Generate problem data
+    rng(1)
+    n = 100;
+    k = 10;
+    F = sprandn(n, k, 0.7);
+    D = sparse(diag( sqrt(k)*rand(n,1) ));
+    mu = randn(n, 1);
+    gamma = 1;
+    Sigma = F*F' + D;
+
+    % Define problem
+    x = sdpvar(n, 1);
+    objective = gamma * (x'*Sigma*x) - mu'*x;
+    constraints = [sum(x) == 1, x >= 0];
+
+    % Solve with OSQP
+    options = sdpsettings('solver', 'osqp');
+    optimize(constraints, objective, options);
+