diff --git a/inst/@sym/simplify.m b/inst/@sym/simplify.m
index 8ce71b2fb..9df5fa698 100644
--- a/inst/@sym/simplify.m
+++ b/inst/@sym/simplify.m
@@ -19,8 +19,13 @@
 %% -*- texinfo -*-
 %% @documentencoding UTF-8
 %% @defmethod @@sym simplify (@var{x})
+%% @defmethodx @@sym simplify (@var{x}, @var{y})
 %% Simplify an expression.
 %%
+%% @var{y} can be 'special' or 'complete'.
+%% 'special' will simplify special function with identities.
+%% 'complete' will replace the identities and then simplify.
+%%
 %% Example:
 %% @example
 %% @group
@@ -47,15 +52,46 @@
 %% @end group
 %% @end example
 %%
+%% When @var{y} is 'special':
+%% @example
+%% @group
+%% syms x
+%% simplify(beta(x, 1), 'special')
+%%   @result{} ans = (sym)
+%%    Γ(x)  
+%%  ────────
+%%  Γ(x + 1)
+%% @end group
+%% @end example
+%%
+%% When @var{y} is 'complete':
+%% @example
+%% @group
+%% syms x
+%% simplify(beta(x, 1), 'complete')
+%%   @result{} ans = (sym)
+%%  1
+%%  ─
+%%  x
+%% @end group
+%% @end example
+%%
 %% @seealso{@@sym/isAlways, @@sym/factor, @@sym/expand, @@sym/rewrite}
 %% @end defmethod
 
+%% Source: http://docs.sympy.org/dev/tutorial/simplification.html
 
-function y = simplify(x)
-
-  cmd = 'return sp.simplify(*_ins),';
+function y = simplify(x, y)
 
-  y = python_cmd (cmd, x);
+  if (nargin == 1)
+    y = python_cmd('return simplify(*_ins),', x);
+  elseif ( strcmp( y, 'special'))
+    y = python_cmd('return expand_func(*_ins),', x);
+  elseif ( strcmp( y, 'complete'))
+    y = python_cmd('return simplify(expand_func(*_ins)),', x);
+  else
+    print_usage ();
+  end
 
 end
 
@@ -66,3 +102,15 @@
 %! q = horner (p);
 %!assert(~isequal( p - q, 0))
 %!assert(isequal( simplify(p - q), 0))
+
+%!test
+%! syms x
+%! A = simplify(beta(x, 1), 'special');
+%! B = gamma(x)/gamma(x+1);
+%! assert( isequal( A, B ))
+
+%!test
+%! syms x
+%! A = simplify(beta(x, 1), 'complete');
+%! B = 1/x;
+%! assert( isequal( A, B ))