diff --git a/advanced/scipy_sparse/examples/direct_solve.py b/advanced/scipy_sparse/examples/direct_solve.py
index e342618bf..38de7cf7e 100644
--- a/advanced/scipy_sparse/examples/direct_solve.py
+++ b/advanced/scipy_sparse/examples/direct_solve.py
@@ -15,7 +15,7 @@
 
 mtx = sp.sparse.lil_array((1000, 1000), dtype=np.float64)
 mtx[0, :100] = rng.random(100)
-mtx[1, 100:200] = mtx[0, :100]
+mtx[1, 100:200] = mtx[[0], :100]
 mtx.setdiag(rng.random(1000))
 
 plt.clf()
@@ -27,4 +27,4 @@
 
 x = sp.sparse.linalg.spsolve(mtx, rhs)
 
-print(f"residual: {np.linalg.norm(mtx * x - rhs)!r}")
+print(f"residual: {np.linalg.norm(mtx @ x - rhs)!r}")
diff --git a/advanced/scipy_sparse/examples/lobpcg_sakurai.py b/advanced/scipy_sparse/examples/lobpcg_sakurai.py
index fc355af75..1ab2618b7 100644
--- a/advanced/scipy_sparse/examples/lobpcg_sakurai.py
+++ b/advanced/scipy_sparse/examples/lobpcg_sakurai.py
@@ -21,7 +21,7 @@ def sakurai(n):
     """
 
     A = sp.sparse.eye(n, n)
-    d0 = np.array(r_[5, 6 * ones(n - 2), 5])
+    d0 = np.hstack([5, 6 * np.ones(n - 2), 5])
     d1 = -4 * np.ones(n)
     d2 = np.ones(n)
     B = sp.sparse.spdiags([d2, d1, d0, d1, d2], [-2, -1, 0, 1, 2], n, n)
@@ -41,13 +41,19 @@ def sakurai(n):
 #
 n = 2500
 A, B, w_ex = sakurai(n)  # Mikota pair
-X = np.rand(n, m)
+X = np.random.random((n, m))
 data = []
-tt = time.clock()
+tt = time.time()
 eigs, vecs, resnh = sp.sparse.linalg.lobpcg(
-    A, X, B, tol=1e-6, maxiter=500, retResidualNormsHistory=1
+    A,
+    X,
+    B,
+    tol=1e-6,
+    largest=False,
+    maxiter=2000,
+    retResidualNormsHistory=1,
 )
-data.append(time.clock() - tt)
+data.append(time.time() - tt)
 print("Results by LOBPCG for n=" + str(n))
 print()
 print(eigs)
diff --git a/advanced/scipy_sparse/examples/pyamg_with_lobpcg.py b/advanced/scipy_sparse/examples/pyamg_with_lobpcg.py
index a7eb4a271..5713a4c29 100644
--- a/advanced/scipy_sparse/examples/pyamg_with_lobpcg.py
+++ b/advanced/scipy_sparse/examples/pyamg_with_lobpcg.py
@@ -7,6 +7,7 @@
 Dirichlet boundary conditions.
 """
 
+import numpy as np
 import scipy as sp
 import matplotlib.pyplot as plt
 
@@ -21,7 +22,7 @@
 ml = smoothed_aggregation_solver(A)
 
 # initial approximation to the K eigenvectors
-X = sp.rand(A.shape[0], K)
+X = np.random.random((A.shape[0], K))
 
 # preconditioner based on ml
 M = ml.aspreconditioner()