diff --git a/.travis.yml b/.travis.yml
index 08425e0db5..5c3e7f50b0 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -15,9 +15,6 @@ matrix:
     - python: 2.7
       env:
         - TEST_SPHINX="true"
-    - python: 2.7
-      env:
-        - TEST_SAGE="true"
     - python: "pypy3"
       env:
         - TEST_DOCTESTS="true"
@@ -27,10 +24,6 @@ matrix:
     - python: "pypy3"
       env:
         - SPLIT="2/2"
-  allow_failures:
-    - python: 2.7
-      env:
-        - TEST_SAGE="true"
 before_install:
   - if [[ "${TEST_SPHINX}" == "true" ]]; then
       sudo apt-get update;
@@ -43,11 +36,6 @@ before_install:
       travis_wait 30 pip install numpy==1.7.1 scipy==0.12.0 gmpy==1.16 Theano==0.6;
       travis_wait 30 pip install --allow-external matplotlib==1.2.1;
     fi
-  - if [[ "${TEST_SAGE}" == "true" ]]; then
-      sudo apt-add-repository -y ppa:aims/sagemath;
-      sudo apt-get update;
-      travis_wait 30 sudo apt-get install sagemath-upstream-binary;
-    fi
   - if [[ "${TEST_PEP}" == "true" ]]; then
       pip install pep8==1.6.2;
     fi
@@ -56,9 +44,7 @@ before_install:
       pip install pytest-cov coveralls --use-mirrors;
     fi
 install:
-  - if [[ "${TEST_SAGE}" != "true" ]]; then
-      python setup.py install;
-    fi
+  python setup.py install;
 script:
   - if [[ "${TEST_PEP}" == "true" ]]; then
       pep8;
diff --git a/bin/test_travis.sh b/bin/test_travis.sh
index 0bf42a2c50..75b1bac092 100755
--- a/bin/test_travis.sh
+++ b/bin/test_travis.sh
@@ -5,9 +5,6 @@ set -e -x # exit on error and echo each command
 if [[ "${TEST_SPHINX}" == "true" ]]; then
     make -C doc html-errors man latex
     LATEXOPTIONS="-interaction=nonstopmode" make -C doc/_build/latex
-elif [[ "${TEST_SAGE}" == "true" ]]; then
-    sage -v
-    sage -python py.test sympy/external/tests/test_sage.py
 else
     if [[ "${TEST_DOCTESTS}" == "true" ]]; then
         cat << EOF | python
diff --git a/sympy/core/add.py b/sympy/core/add.py
index 539225bad1..9a02e620bd 100644
--- a/sympy/core/add.py
+++ b/sympy/core/add.py
@@ -725,12 +725,6 @@ def _eval_transpose(self):
     def __neg__(self):
         return self.func(*[-t for t in self.args])
 
-    def _sage_(self):
-        s = 0
-        for x in self.args:
-            s += x._sage_()
-        return s
-
     def primitive(self):
         """
         Return ``(R, self/R)`` where ``R``` is the Rational GCD of ``self```.
diff --git a/sympy/core/function.py b/sympy/core/function.py
index 89c0e211ec..ed9380d991 100644
--- a/sympy/core/function.py
+++ b/sympy/core/function.py
@@ -676,13 +676,6 @@ def _eval_as_leading_term(self, x):
         else:
             return self.func(*args)
 
-    def _sage_(self):
-        import sage.all as sage
-        fname = self.func.__name__
-        func = getattr(sage, fname)
-        args = [arg._sage_() for arg in self.args]
-        return func(*args)
-
 
 class AppliedUndef(Function):
     """
@@ -698,12 +691,6 @@ def __new__(cls, *args, **options):
     def _eval_as_leading_term(self, x):
         return self
 
-    def _sage_(self):
-        import sage.all as sage
-        fname = str(self.func)
-        args = [arg._sage_() for arg in self.args]
-        func = sage.function(fname, *args)
-        return func
 
 class UndefinedFunction(FunctionClass):
     """
@@ -1320,11 +1307,6 @@ def _eval_nseries(self, x, n, logx):
     def _eval_as_leading_term(self, x):
         return self.args[0].as_leading_term(x)
 
-    def _sage_(self):
-        import sage.all as sage
-        args = [arg._sage_() for arg in self.args]
-        return sage.derivative(*args)
-
 
 class Lambda(Expr):
     """
diff --git a/sympy/core/mul.py b/sympy/core/mul.py
index 89addd65f7..5955b6516e 100644
--- a/sympy/core/mul.py
+++ b/sympy/core/mul.py
@@ -1465,12 +1465,6 @@ def _eval_transpose(self):
     def _eval_adjoint(self):
         return self.func(*[t.adjoint() for t in self.args[::-1]])
 
-    def _sage_(self):
-        s = 1
-        for x in self.args:
-            s *= x._sage_()
-        return s
-
     def as_content_primitive(self, radical=False):
         """Return the tuple (R, self/R) where R is the positive Rational
         extracted from self.
diff --git a/sympy/core/numbers.py b/sympy/core/numbers.py
index c9a0680ed6..c37ae64934 100644
--- a/sympy/core/numbers.py
+++ b/sympy/core/numbers.py
@@ -1032,10 +1032,6 @@ def __hash__(self):
     def epsilon_eq(self, other, epsilon="1e-15"):
         return abs(self - other) < Float(epsilon)
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.RealNumber(str(self))
-
     def __format__(self, format_spec):
         return format(decimal.Decimal(str(self)), format_spec)
 
@@ -1557,10 +1553,6 @@ def lcm(self, other):
     def as_numer_denom(self):
         return Integer(self.p), Integer(self.q)
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.Integer(self.p)/sage.Integer(self.q)
-
     def as_content_primitive(self, radical=False):
         """Return the tuple (R, self/R) where R is the positive Rational
         extracted from self.
@@ -2471,10 +2463,6 @@ def _eval_power(self, expt):
     def _as_mpf_val(self, prec):
         return mlib.finf
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.oo
-
     def __hash__(self):
         return super(Infinity, self).__hash__()
 
@@ -2682,10 +2670,6 @@ def _eval_power(self, expt):
     def _as_mpf_val(self, prec):
         return mlib.fninf
 
-    def _sage_(self):
-        import sage.all as sage
-        return -(sage.oo)
-
     def __hash__(self):
         return super(NegativeInfinity, self).__hash__()
 
@@ -2835,10 +2819,6 @@ def __div__(self, other):
     def _as_mpf_val(self, prec):
         return _mpf_nan
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.NaN
-
     def __hash__(self):
         return super(NaN, self).__hash__()
 
@@ -2927,10 +2907,6 @@ def _eval_power(self, expt):
                 else:
                     return S.Zero
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.UnsignedInfinityRing.gen()
-
 
 zoo = S.ComplexInfinity
 
@@ -3107,9 +3083,6 @@ def _eval_rewrite_as_cos(self):
         I = S.ImaginaryUnit
         return cos(I) + I*cos(I + S.Pi/2)
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.e
 E = S.Exp1
 
 
@@ -3175,9 +3148,6 @@ def approximation_interval(self, number_cls):
         elif issubclass(number_cls, Rational):
             return (Rational(223, 71), Rational(22, 7))
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.pi
 pi = S.Pi
 
 
@@ -3238,10 +3208,6 @@ def approximation_interval(self, number_cls):
         elif issubclass(number_cls, Rational):
             pass
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.golden_ratio
-
 
 class EulerGamma(with_metaclass(Singleton, NumberSymbol)):
     r"""The Euler-Mascheroni constant.
@@ -3298,10 +3264,6 @@ def approximation_interval(self, number_cls):
         elif issubclass(number_cls, Rational):
             return (S.Half, Rational(3, 5))
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.euler_gamma
-
 
 class Catalan(with_metaclass(Singleton, NumberSymbol)):
     r"""Catalan's constant.
@@ -3351,10 +3313,6 @@ def approximation_interval(self, number_cls):
         elif issubclass(number_cls, Rational):
             return (Rational(9, 10), S.One)
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.catalan
-
 
 class ImaginaryUnit(with_metaclass(Singleton, AtomicExpr)):
     r"""The imaginary unit, `i = \sqrt{-1}`.
@@ -3431,10 +3389,6 @@ def _eval_power(self, expt):
     def as_base_exp(self):
         return S.NegativeOne, S.Half
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.I
-
 I = S.ImaginaryUnit
 
 
diff --git a/sympy/core/power.py b/sympy/core/power.py
index 0655974f57..31f50ac379 100644
--- a/sympy/core/power.py
+++ b/sympy/core/power.py
@@ -1317,9 +1317,6 @@ def _taylor_term(self, n, x, *previous_terms): # of (1+x)**e
         from sympy import binomial
         return binomial(self.exp, n) * self.func(x, n)
 
-    def _sage_(self):
-        return self.args[0]._sage_()**self.args[1]._sage_()
-
     def as_content_primitive(self, radical=False):
         """Return the tuple (R, self/R) where R is the positive Rational
         extracted from self.
diff --git a/sympy/core/symbol.py b/sympy/core/symbol.py
index eb8d802e1c..8bd4173858 100644
--- a/sympy/core/symbol.py
+++ b/sympy/core/symbol.py
@@ -148,10 +148,6 @@ def as_real_imag(self, deep=True, **hints):
         else:
             return (re(self), im(self))
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.var(self.name)
-
     def is_constant(self, *wrt, **flags):
         if not wrt:
             return False
diff --git a/sympy/core/tests/test_sympify.py b/sympy/core/tests/test_sympify.py
index af7b3ced92..ba52f4ea82 100644
--- a/sympy/core/tests/test_sympify.py
+++ b/sympy/core/tests/test_sympify.py
@@ -212,19 +212,6 @@ def test_sympify_factorial():
     raises(SympifyError, lambda: sympify("x!!!"))
 
 
-def test_sage():
-    # how to effectivelly test for the _sage_() method without having SAGE
-    # installed?
-    assert hasattr(x, "_sage_")
-    assert hasattr(Integer(3), "_sage_")
-    assert hasattr(sin(x), "_sage_")
-    assert hasattr(cos(x), "_sage_")
-    assert hasattr(x**2, "_sage_")
-    assert hasattr(x + y, "_sage_")
-    assert hasattr(exp(x), "_sage_")
-    assert hasattr(log(x), "_sage_")
-
-
 def test_issue_3595():
     assert sympify("a_") == Symbol("a_")
     assert sympify("_a") == Symbol("_a")
diff --git a/sympy/external/tests/test_sage.py b/sympy/external/tests/test_sage.py
deleted file mode 100644
index d2dcfab5ae..0000000000
--- a/sympy/external/tests/test_sage.py
+++ /dev/null
@@ -1,254 +0,0 @@
-# This testfile tests SymPy <-> Sage compatibility
-#
-# Execute this test inside Sage, e.g. with:
-# sage -python bin/test sympy/external/tests/test_sage.py
-#
-# This file can be tested by Sage itself by:
-# sage -t sympy/external/tests/test_sage.py
-# and if all tests pass, it should be copied (verbatim) to Sage, so that it is
-# automatically doctested by Sage.  Note that this second method imports the
-# version of SymPy in Sage, whereas the -python method imports the local version
-# of SymPy (both use the local version of the tests, however).
-#
-# Don't test any SymPy features here. Just pure interaction with Sage.
-# Always write regular SymPy tests for anything, that can be tested in pure
-# Python (without Sage). Here we test everything, that a user may need when
-# using SymPy with Sage.
-
-import os
-import re
-import sys
-
-from sympy.external import import_module
-
-sage = import_module('sage.all', __import__kwargs={'fromlist': ['all']})
-if not sage:
-    #bin/test will not execute any tests now
-    disabled = True
-
-if sys.version_info[0] == 3:
-    # Sage does not support Python 3 currently
-    disabled = True
-
-import sympy
-
-from sympy.utilities.pytest import XFAIL
-
-
-def check_expression(expr, var_symbols, only_from_sympy=False):
-    """
-    Does eval(expr) both in Sage and SymPy and does other checks.
-    """
-
-    # evaluate the expression in the context of Sage:
-    if var_symbols:
-        sage.var(var_symbols)
-    a = globals().copy()
-    # safety checks...
-    a.update(sage.__dict__)
-    assert "sin" in a
-    is_different = False
-    try:
-        e_sage = eval(expr, a)
-        assert not isinstance(e_sage, sympy.Basic)
-    except (NameError, TypeError):
-        is_different = True
-        pass
-
-    # evaluate the expression in the context of SymPy:
-    if var_symbols:
-        sympy_vars = sympy.var(var_symbols)
-    b = globals().copy()
-    b.update(sympy.__dict__)
-    assert "sin" in b
-    b.update(sympy.__dict__)
-    e_sympy = eval(expr, b)
-    assert isinstance(e_sympy, sympy.Basic)
-
-    # Sympy func may have specific _sage_ method
-    if is_different:
-        _sage_method = getattr(e_sympy.func, "_sage_")
-        e_sage = _sage_method(sympy.S(e_sympy))
-
-    # Do the actual checks:
-    if not only_from_sympy:
-        assert sympy.S(e_sage) == e_sympy
-    assert e_sage == sage.SR(e_sympy)
-
-
-def test_basics():
-    check_expression("x", "x")
-    check_expression("x**2", "x")
-    check_expression("x**2+y**3", "x y")
-    check_expression("1/(x+y)**2-x**3/4", "x y")
-
-
-def test_complex():
-    check_expression("I", "")
-    check_expression("23+I*4", "x")
-
-
-@XFAIL
-def test_complex_fail():
-    # Sage doesn't properly implement _sympy_ on I
-    check_expression("I*y", "y")
-    check_expression("x+I*y", "x y")
-
-
-def test_integer():
-    check_expression("4*x", "x")
-    check_expression("-4*x", "x")
-
-
-def test_real():
-    check_expression("1.123*x", "x")
-    check_expression("-18.22*x", "x")
-
-
-def test_E():
-    assert sympy.sympify(sage.e) == sympy.E
-    assert sage.e == sage.SR(sympy.E)
-
-
-def test_pi():
-    assert sympy.sympify(sage.pi) == sympy.pi
-    assert sage.pi == sage.SR(sympy.pi)
-
-
-def test_euler_gamma():
-    assert sympy.sympify(sage.euler_gamma) == sympy.EulerGamma
-    assert sage.euler_gamma == sage.SR(sympy.EulerGamma)
-
-
-def test_oo():
-    assert sympy.sympify(sage.oo) == sympy.oo
-    assert sage.oo == sage.SR(sympy.oo)
-    assert sympy.sympify(-sage.oo) == -sympy.oo
-    assert -sage.oo == sage.SR(-sympy.oo)
-    #assert sympy.sympify(sage.UnsignedInfinityRing.gen()) == sympy.zoo
-    #assert sage.UnsignedInfinityRing.gen() == sage.SR(sympy.zoo)
-
-def test_NaN():
-    assert sympy.sympify(sage.NaN) == sympy.nan
-    assert sage.NaN == sage.SR(sympy.nan)
-
-
-def test_Catalan():
-    assert sympy.sympify(sage.catalan) == sympy.Catalan
-    assert sage.catalan == sage.SR(sympy.Catalan)
-
-
-def test_GoldenRation():
-    assert sympy.sympify(sage.golden_ratio) == sympy.GoldenRatio
-    assert sage.golden_ratio == sage.SR(sympy.GoldenRatio)
-
-
-def test_functions():
-    # Test at least one Function without own _sage_ method
-    assert not "_sage_" in sympy.factorial.__dict__
-    check_expression("factorial(x)", "x")
-    check_expression("sin(x)", "x")
-    check_expression("cos(x)", "x")
-    check_expression("tan(x)", "x")
-    check_expression("cot(x)", "x")
-    check_expression("asin(x)", "x")
-    check_expression("acos(x)", "x")
-    check_expression("atan(x)", "x")
-    check_expression("atan2(y, x)", "x, y")
-    check_expression("acot(x)", "x")
-    check_expression("sinh(x)", "x")
-    check_expression("cosh(x)", "x")
-    check_expression("tanh(x)", "x")
-    check_expression("coth(x)", "x")
-    check_expression("asinh(x)", "x")
-    check_expression("acosh(x)", "x")
-    check_expression("atanh(x)", "x")
-    check_expression("acoth(x)", "x")
-    check_expression("exp(x)", "x")
-    check_expression("log(x)", "x")
-    check_expression("re(x)", "x")
-    check_expression("im(x)", "x")
-    check_expression("sign(x)", "x")
-    check_expression("abs(x)", "x")
-    check_expression("arg(x)", "x")
-    check_expression("conjugate(x)", "x")
-
-    # The following tests differently named functions
-    check_expression("besselj(y, x)", "x, y")
-    check_expression("bessely(y, x)", "x, y")
-    check_expression("besseli(y, x)", "x, y")
-    check_expression("besselk(y, x)", "x, y")
-    check_expression("DiracDelta(x)", "x")
-    check_expression("KroneckerDelta(x, y)", "x, y")
-    check_expression("expint(y, x)", "x, y")
-    check_expression("Si(x)", "x")
-    check_expression("Ci(x)", "x")
-    check_expression("Shi(x)", "x")
-    check_expression("Chi(x)", "x")
-    check_expression("loggamma(x)", "x")
-    check_expression("Ynm(n,m,x,y)", "n, m, x, y")
-    check_expression("hyper((n,m),(m,n),x)", "n, m, x")
-
-def test_issue_4023():
-    sage.var("a x")
-    log = sage.log
-    i = sympy.integrate(log(x)/a, (x, a, a + 1))
-    i2 = sympy.simplify(i)
-    s = sage.SR(i2)
-    assert s == (a*log(1 + a) - a*log(a) + log(1 + a) - 1)/a
-
-def test_integral():
-    #test Sympy-->Sage
-    check_expression("Integral(x, (x,))", "x", only_from_sympy=True)
-    check_expression("Integral(x, (x, 0, 1))", "x", only_from_sympy=True)
-    check_expression("Integral(x*y, (x,), (y, ))", "x,y", only_from_sympy=True)
-    check_expression("Integral(x*y, (x,), (y, 0, 1))", "x,y", only_from_sympy=True)
-    check_expression("Integral(x*y, (x, 0, 1), (y,))", "x,y", only_from_sympy=True)
-    check_expression("Integral(x*y, (x, 0, 1), (y, 0, 1))", "x,y", only_from_sympy=True)
-    check_expression("Integral(x*y*z, (x, 0, 1), (y, 0, 1), (z, 0, 1))", "x,y,z", only_from_sympy=True)
-
-@XFAIL
-def test_integral_failing():
-    check_expression("Integral(x, (x, 0))", "x", only_from_sympy=True)
-    check_expression("Integral(x*y, (x,), (y, 0))", "x,y", only_from_sympy=True)
-    check_expression("Integral(x*y, (x, 0, 1), (y, 0))", "x,y", only_from_sympy=True)
-
-def test_undefined_function():
-    f = sympy.Function('f')
-    sf = sage.function('f')
-    x = sympy.symbols('x')
-    sx = sage.var('x')
-    assert bool(sf(sx) == f(x)._sage_())
-    #assert bool(f == sympy.sympify(sf))
-
-# This string contains Sage doctests, that execute all the functions above.
-# When you add a new function, please add it here as well.
-"""
-
-TESTS::
-
-    sage: test_basics()
-    sage: test_basics()
-    sage: test_complex()
-    sage: test_integer()
-    sage: test_real()
-    sage: test_E()
-    sage: test_pi()
-    sage: test_euler_gamma()
-    sage: test_oo()
-    sage: test_NaN()
-    sage: test_Catalan()
-    sage: test_GoldenRation()
-    sage: test_functions()
-    sage: test_issue_4023()
-    sage: test_integral()
-    sage: test_undefined_function()
-
-Sage has no symbolic Lucas function at the moment::
-
-    sage: check_expression("lucas(x)", "x")
-    Traceback (most recent call last):
-    ...
-    AttributeError: 'module' object has no attribute 'lucas'
-
-"""
diff --git a/sympy/functions/combinatorial/factorials.py b/sympy/functions/combinatorial/factorials.py
index 0ff92c6fb2..399d7c9456 100644
--- a/sympy/functions/combinatorial/factorials.py
+++ b/sympy/functions/combinatorial/factorials.py
@@ -470,10 +470,6 @@ def _eval_is_integer(self):
         return fuzzy_and((self.args[0].is_integer, self.args[1].is_integer,
                           self.args[1].is_nonnegative))
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.rising_factorial(self.args[0]._sage_(), self.args[1]._sage_())
-
 
 class FallingFactorial(CombinatorialFunction):
     """Falling factorial (related to rising factorial) is a double valued
@@ -545,10 +541,6 @@ def _eval_is_integer(self):
         return fuzzy_and((self.args[0].is_integer, self.args[1].is_integer,
                           self.args[1].is_nonnegative))
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.falling_factorial(self.args[0]._sage_(), self.args[1]._sage_())
-
 
 rf = RisingFactorial
 ff = FallingFactorial
diff --git a/sympy/functions/elementary/complexes.py b/sympy/functions/elementary/complexes.py
index 87e68f5c68..8fbd619d76 100644
--- a/sympy/functions/elementary/complexes.py
+++ b/sympy/functions/elementary/complexes.py
@@ -105,10 +105,6 @@ def _eval_rewrite_as_im(self, arg):
     def _eval_is_algebraic(self):
         return self.args[0].is_algebraic
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.real_part(self.args[0]._sage_())
-
 
 class im(Function):
     """
@@ -202,10 +198,6 @@ def _eval_derivative(self, x):
             return -S.ImaginaryUnit \
                 * re(Derivative(self.args[0], x, evaluate=True))
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.imag_part(self.args[0]._sage_())
-
     def _eval_rewrite_as_re(self, arg):
         return self.args[0] - re(self.args[0])
 
@@ -357,10 +349,6 @@ def _eval_power(self, other):
         ):
             return S.One
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.sgn(self.args[0]._sage_())
-
     def _eval_rewrite_as_Piecewise(self, arg):
         if arg.is_extended_real:
             return Piecewise((1, arg > 0), (-1, arg < 0), (0, True))
@@ -539,10 +527,6 @@ def _eval_nseries(self, x, n, logx):
             (sign(direction)*s, True),
         )
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.abs_symbolic(self.args[0]._sage_())
-
     def _eval_derivative(self, x):
         if self.args[0].is_extended_real or self.args[0].is_imaginary:
             return Derivative(self.args[0], x, evaluate=True) \
diff --git a/sympy/functions/elementary/hyperbolic.py b/sympy/functions/elementary/hyperbolic.py
index 6d7264528a..6eb03560e3 100644
--- a/sympy/functions/elementary/hyperbolic.py
+++ b/sympy/functions/elementary/hyperbolic.py
@@ -736,10 +736,6 @@ def taylor_term(n, x, *previous_terms):
     def _eval_rewrite_as_cosh(self, arg):
         return S.ImaginaryUnit / cosh(arg + S.ImaginaryUnit * S.Pi / 2)
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.csch(self.args[0]._sage_())
-
 
 class sech(ReciprocalHyperbolicFunction):
     r"""
@@ -781,11 +777,6 @@ def taylor_term(n, x, *previous_terms):
     def _eval_rewrite_as_sinh(self, arg):
         return S.ImaginaryUnit / sinh(arg + S.ImaginaryUnit * S.Pi /2)
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.sech(self.args[0]._sage_())
-
-
 
 ###############################################################################
 ############################# HYPERBOLIC INVERSES #############################
diff --git a/sympy/functions/special/bessel.py b/sympy/functions/special/bessel.py
index 50f19e23da..38b966f2aa 100644
--- a/sympy/functions/special/bessel.py
+++ b/sympy/functions/special/bessel.py
@@ -197,10 +197,6 @@ def _eval_is_extended_real(self):
         if nu.is_integer and z.is_extended_real:
             return True
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.bessel_J(self.args[0]._sage_(), self.args[1]._sage_())
-
 
 class bessely(BesselBase):
     r"""
@@ -276,10 +272,6 @@ def _eval_is_extended_real(self):
         if nu.is_integer and z.is_positive:
             return True
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.bessel_Y(self.args[0]._sage_(), self.args[1]._sage_())
-
 
 class besseli(BesselBase):
     r"""
@@ -376,10 +368,6 @@ def _eval_is_extended_real(self):
         if nu.is_integer and z.is_extended_real:
             return True
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.bessel_I(self.args[0]._sage_(), self.args[1]._sage_())
-
 
 class besselk(BesselBase):
     r"""
@@ -460,10 +448,6 @@ def _eval_is_extended_real(self):
         if nu.is_integer and z.is_positive:
             return True
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.bessel_K(self.args[0]._sage_(), self.args[1]._sage_())
-
 
 class hankel1(BesselBase):
     r"""
diff --git a/sympy/functions/special/delta_functions.py b/sympy/functions/special/delta_functions.py
index 4dedaa20a3..7ccebe0d98 100644
--- a/sympy/functions/special/delta_functions.py
+++ b/sympy/functions/special/delta_functions.py
@@ -164,10 +164,6 @@ def is_simple(self, x):
     def _latex_no_arg(printer):
         return r'\delta'
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.dirac_delta(self.args[0]._sage_())
-
 
 ###############################################################################
 ############################## HEAVISIDE FUNCTION #############################
@@ -250,7 +246,3 @@ def _eval_rewrite_as_Piecewise(self, arg):
     def _eval_rewrite_as_sign(self, arg):
         if arg.is_extended_real:
             return (sign(arg)+1)/2
-
-    def _sage_(self):
-        import sage.all as sage
-        return sage.heaviside(self.args[0]._sage_())
diff --git a/sympy/functions/special/elliptic_integrals.py b/sympy/functions/special/elliptic_integrals.py
index 53131139fb..9917317808 100644
--- a/sympy/functions/special/elliptic_integrals.py
+++ b/sympy/functions/special/elliptic_integrals.py
@@ -80,10 +80,6 @@ def _eval_rewrite_as_hyper(self, z):
     def _eval_rewrite_as_meijerg(self, z):
         return meijerg(((S.Half, S.Half), []), ((S.Zero,), (S.Zero,)), -z)/2
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.elliptic_kc(self.args[0]._sage_())
-
 
 class elliptic_f(Function):
     r"""
diff --git a/sympy/functions/special/error_functions.py b/sympy/functions/special/error_functions.py
index 1beededf15..5d4fb294e7 100644
--- a/sympy/functions/special/error_functions.py
+++ b/sympy/functions/special/error_functions.py
@@ -1277,10 +1277,6 @@ def _eval_nseries(self, x, n, logx):
                 return f._eval_nseries(x, n, logx)
         return super(expint, self)._eval_nseries(x, n, logx)
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.exp_integral_e(self.args[0]._sage_(), self.args[1]._sage_())
-
 
 def E1(z):
     """
@@ -1682,9 +1678,6 @@ def _eval_rewrite_as_expint(self, z):
         # XXX should we polarify z?
         return pi/2 + (E1(polar_lift(I)*z) - E1(polar_lift(-I)*z))/2/I
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.sin_integral(self.args[0]._sage_())
 
 class Ci(TrigonometricIntegral):
     r"""
@@ -1784,10 +1777,6 @@ def _Ifactor(cls, z, sign):
     def _eval_rewrite_as_expint(self, z):
         return -(E1(polar_lift(I)*z) + E1(polar_lift(-I)*z))/2
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.cos_integral(self.args[0]._sage_())
-
 
 class Shi(TrigonometricIntegral):
     r"""
@@ -1873,10 +1862,6 @@ def _eval_rewrite_as_expint(self, z):
         # XXX should we polarify z?
         return (E1(z) - E1(exp_polar(I*pi)*z))/2 - I*pi/2
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.sinh_integral(self.args[0]._sage_())
-
 
 class Chi(TrigonometricIntegral):
     r"""
@@ -1986,10 +1971,6 @@ def _latex(self, printer, exp=None):
     def _latex_no_arg(printer):
         return r'\operatorname{Chi}'
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.cosh_integral(self.args[0]._sage_())
-
 
 ###############################################################################
 #################### FRESNEL INTEGRALS ########################################
diff --git a/sympy/functions/special/gamma_functions.py b/sympy/functions/special/gamma_functions.py
index 719617ae4c..4735c84fcc 100644
--- a/sympy/functions/special/gamma_functions.py
+++ b/sympy/functions/special/gamma_functions.py
@@ -936,10 +936,6 @@ def fdiff(self, argindex=1):
         else:
             raise ArgumentIndexError(self, argindex)
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.log_gamma(self.args[0]._sage_())
-
 
 def digamma(x):
     r"""
diff --git a/sympy/functions/special/hyper.py b/sympy/functions/special/hyper.py
index 94a1f2ea6d..2d57882f2f 100644
--- a/sympy/functions/special/hyper.py
+++ b/sympy/functions/special/hyper.py
@@ -305,12 +305,6 @@ def _eval_simplify(self, ratio, measure):
         from sympy.simplify.hyperexpand import hyperexpand
         return hyperexpand(self)
 
-    def _sage_(self):
-        import sage.all as sage
-        ap = [arg._sage_() for arg in self.args[0]]
-        bq = [arg._sage_() for arg in self.args[1]]
-        return sage.hypergeometric(ap, bq, self.argument._sage_())
-
 
 class meijerg(TupleParametersBase):
     r"""
diff --git a/sympy/functions/special/spherical_harmonics.py b/sympy/functions/special/spherical_harmonics.py
index 7d905b6110..5db886678e 100644
--- a/sympy/functions/special/spherical_harmonics.py
+++ b/sympy/functions/special/spherical_harmonics.py
@@ -231,13 +231,6 @@ def _eval_evalf(self, prec):
             res = mp.spherharm(n, m, theta, phi)
         return Expr._from_mpmath(res, prec)
 
-    def _sage_(self):
-        import sage.all as sage
-        return sage.spherical_harmonic(self.args[0]._sage_(),
-                                       self.args[1]._sage_(),
-                                       self.args[2]._sage_(),
-                                       self.args[3]._sage_())
-
 
 def Ynm_c(n, m, theta, phi):
     r"""Conjugate spherical harmonics defined as
diff --git a/sympy/functions/special/tensor_functions.py b/sympy/functions/special/tensor_functions.py
index 82512cf61b..53d7a4c685 100644
--- a/sympy/functions/special/tensor_functions.py
+++ b/sympy/functions/special/tensor_functions.py
@@ -445,7 +445,3 @@ def _get_preferred_index(self):
     @staticmethod
     def _latex_no_arg(printer):
         return r'\delta'
-
-    def _sage_(self):
-        import sage.all as sage
-        return sage.kronecker_delta(self.args[0]._sage_(), self.args[1]._sage_())
diff --git a/sympy/integrals/integrals.py b/sympy/integrals/integrals.py
index dbf05c9a3e..5ee8bd6a5a 100644
--- a/sympy/integrals/integrals.py
+++ b/sympy/integrals/integrals.py
@@ -1067,30 +1067,6 @@ def as_sum(self, n, method="midpoint"):
             result += self.function.subs(sym, xi)
         return result*dx
 
-    def _sage_(self):
-        import sage.all as sage
-        f, limits = self.function._sage_(), list(self.limits)
-        for limit in limits:
-            if len(limit) == 1:
-                x = limit[0]
-                f = sage.integral(f,
-                                    x._sage_(),
-                                    hold=True)
-            elif len(limit) == 2:
-                x, b = limit
-                f = sage.integral(f,
-                                    x._sage_(),
-                                    b._sage_(),
-                                    hold=True)
-            else:
-                x, a, b = limit
-                f = sage.integral(f,
-                                  (x._sage_(),
-                                    a._sage_(),
-                                    b._sage_()),
-                                    hold=True)
-        return f
-
 
 @xthreaded
 def integrate(*args, **kwargs):
diff --git a/sympy/series/order.py b/sympy/series/order.py
index 343c407901..c55b345b41 100644
--- a/sympy/series/order.py
+++ b/sympy/series/order.py
@@ -413,8 +413,4 @@ def _eval_transpose(self):
         if expr is not None:
             return self.func(expr, *self.args[1:])
 
-    def _sage_(self):
-        #XXX: SAGE doesn't have Order yet. Let's return 0 instead.
-        return Rational(0)._sage_()
-
 O = Order
diff --git a/sympy/utilities/runtests.py b/sympy/utilities/runtests.py
index f0fbc98ed1..325ade52ac 100644
--- a/sympy/utilities/runtests.py
+++ b/sympy/utilities/runtests.py
@@ -223,8 +223,7 @@ def run_all_tests(test_args=(), test_kwargs={}, doctest_args=(),
     Run all tests.
 
     Right now, this runs the regular tests (py.test), the doctests
-    (bin/doctest), the examples (examples/all.py), and the sage tests (see
-    sympy/external/tests/test_sage.py).
+    (bin/doctest), the examples (examples/all.py).
 
     This is what ``setup.py test`` uses.
 
@@ -258,16 +257,6 @@ def run_all_tests(test_args=(), test_kwargs={}, doctest_args=(),
         if not run_examples(*examples_args, **examples_kwargs):
             tests_successful = False
 
-        # Sage tests
-        if not (sys.platform == "win32" or PY3):
-            # run Sage tests; Sage currently doesn't support Windows or Python 3
-            dev_null = open(os.devnull, 'w')
-            if subprocess.call("sage -v", shell=True, stdout=dev_null,
-                               stderr=dev_null) == 0:
-                if subprocess.call("sage -python py.test "
-                                   "sympy/external/tests/test_sage.py", shell=True) != 0:
-                    tests_successful = False
-
         if tests_successful:
             return
         else: