diff --git a/Violin.m b/Violin.m
index e1c74ca..bd3c9a3 100644
--- a/Violin.m
+++ b/Violin.m
@@ -106,13 +106,13 @@
             args = obj.checkInputs(data, pos, varargin{:});
             data = data(not(isnan(data)));
             if numel(data) == 1
-                obj.MedianPlot = scatter(pos, data, 'filled');
+                obj.MedianPlot = scatter(args.Parent, pos, data, 'filled');
                 obj.MedianColor = args.MedianColor;
                 obj.MedianPlot.MarkerEdgeColor = args.EdgeColor;
                 return
             end
 
-            hold('on');
+            hold(args.Parent, 'on');
 
             % calculate kernel density estimation for the violin
             if isempty(data)
@@ -139,19 +139,19 @@
             end
             jitter = 2*(rand(size(data))-0.5);
             obj.ScatterPlot = ...
-                scatter(pos + jitter.*jitterstrength, data, 'filled');
+                scatter(args.Parent, pos + jitter.*jitterstrength, data, 'filled');
 
             % plot the violin
             obj.ViolinPlot =  ... % plot color will be overwritten later
                 fill([pos+density*width pos-density(end:-1:1)*width], ...
-                     [value value(end:-1:1)], [1 1 1]);
+                     [value value(end:-1:1)], [1 1 1], 'Parent', args.Parent);
 
             % plot the mini-boxplot within the violin
             quartiles = quantile(data, [0.25, 0.5, 0.75]);         
             obj.BoxPlot = ... % plot color will be overwritten later
                 fill(pos+[-1,1,1,-1]*args.BoxWidth, ...
                      [quartiles(1) quartiles(1) quartiles(3) quartiles(3)], ...
-                     [1 1 1]);
+                     [1 1 1], 'Parent', args.Parent);
                  
             % plot the data mean
             meanValue = mean(data);
@@ -163,7 +163,7 @@
             if meanDensityWidth<args.BoxWidth/2
                 meanDensityWidth=args.BoxWidth/2;
             end
-            obj.MeanPlot = plot(pos+[-1,1].*meanDensityWidth, ...
+            obj.MeanPlot = plot(args.Parent, pos+[-1,1].*meanDensityWidth, ...
                                 [meanValue, meanValue]);
             obj.MeanPlot.LineWidth = 1;
                  
@@ -173,15 +173,15 @@
             hiwhisker = quartiles(3) + 1.5*IQR;
             hiwhisker = min(hiwhisker, max(data(data < hiwhisker)));
             if ~isempty(lowhisker) && ~isempty(hiwhisker)
-                obj.WhiskerPlot = plot([pos pos], [lowhisker hiwhisker]);
+                obj.WhiskerPlot = plot(args.Parent, [pos pos], [lowhisker hiwhisker]);
             end
-            obj.MedianPlot = scatter(pos, quartiles(2), [], [1 1 1], 'filled');
+            obj.MedianPlot = scatter(args.Parent, pos, quartiles(2), [], [1 1 1], 'filled');
 
             obj.NotchPlots = ...
-                 scatter(pos, quartiles(2)-1.57*IQR/sqrt(length(data)), ...
+                 scatter(args.Parent, pos, quartiles(2)-1.57*IQR/sqrt(length(data)), ...
                          [], [1 1 1], 'filled', '^');
             obj.NotchPlots(2) = ...
-                 scatter(pos, quartiles(2)+1.57*IQR/sqrt(length(data)), ...
+                 scatter(args.Parent, pos, quartiles(2)+1.57*IQR/sqrt(length(data)), ...
                          [], [1 1 1], 'filled', 'v');
 
             obj.EdgeColor = args.EdgeColor;
@@ -197,6 +197,7 @@
             obj.ShowData = args.ShowData;
             obj.ShowNotches = args.ShowNotches;
             obj.ShowMean = args.ShowMean;
+            uistack(obj.ViolinPlot,'bottom');
         end
 
         function set.EdgeColor(obj, color)
@@ -338,6 +339,7 @@
             p.addParameter('ShowData', true, isscalarlogical);
             p.addParameter('ShowNotches', false, isscalarlogical);
             p.addParameter('ShowMean', false, isscalarlogical);
+            p.addParameter('Parent', gca, @ishghandle);
 
             p.parse(data, pos, varargin{:});
             results = p.Results;
diff --git a/violinplot.m b/violinplot.m
index 57264f8..c5f5c96 100644
--- a/violinplot.m
+++ b/violinplot.m
@@ -1,4 +1,4 @@
-function violins = violinplot(data, cats, varargin)
+function varargout = violinplot(data, cats, varargin)
 %Violinplots plots violin plots of some data and categories
 %   VIOLINPLOT(DATA) plots a violin of a double vector DATA
 %
@@ -46,17 +46,19 @@
 %                    Defaults to false
 %     'GroupOrder'   Cell of category names in order to be plotted.
 %                    Defaults to alphabetical ordering
+%     'Parent'       Axis handle to plot in
+%                    Defaults to gca
 
 % Copyright (c) 2016, Bastian Bechtold
 % This code is released under the terms of the BSD 3-clause license
 
     hascategories = exist('cats','var') && not(isempty(cats));
     
-    %parse the optional grouporder argument 
-    %if it exists parse the categories order 
+    % parse the optional grouporder argument 
+    % if it exists parse the categories order 
     % but also delete it from the arguments passed to Violin
     grouporder = {};
-    idx=find(strcmp(varargin, 'GroupOrder'));
+    idx = find(strcmp(varargin, 'GroupOrder'));
     if ~isempty(idx) && numel(varargin)>idx
         if iscell(varargin{idx+1})
             grouporder = varargin{idx+1};
@@ -65,6 +67,19 @@
             error('Second argument of ''GroupOrder'' optional arg must be a cell of category names')
         end
     end
+    
+    % parse the optional axis handle
+    AxisHandle = gca;
+    idx = find(strcmp(varargin, 'Parent'));
+    if ~isempty(idx) && numel(varargin)>idx
+        if ishghandle(varargin{idx+1})
+            AxisHandle = varargin{idx+1};
+        else
+            error('Second argument of ''Parent'' optional arg must be a valid axis handle')
+        end
+    end
+    
+    
 
     % tabular data
     if isa(data, 'dataset') || isstruct(data) || istable(data)
@@ -85,7 +100,7 @@
             thisData = data.(catnames{n});
             violins(n) = Violin(thisData, n, varargin{:});
         end
-        set(gca, 'XTick', 1:length(catnames), 'XTickLabels', catnames);
+        set(AxisHandle, 'XTick', 1:length(catnames), 'XTickLabel', catnames);
 
     % 1D data, one category for each data point
     elseif hascategories && numel(data) == numel(cats)
@@ -104,12 +119,12 @@
             thisData = data(cats == thisCat);
             violins(n) = Violin(thisData, n, varargin{:});
         end
-        set(gca, 'XTick', 1:length(catnames), 'XTickLabels', catnames_labels);
+        set(AxisHandle, 'XTick', 1:length(catnames), 'XTickLabel', catnames_labels);
 
     % 1D data, no categories
     elseif not(hascategories) && isvector(data)
         violins = Violin(data, 1, varargin{:});
-        set(gca, 'XTick', 1);
+        set(AxisHandle, 'XTick', 1);
 
     % 2D data with or without categories
     elseif ismatrix(data)
@@ -117,11 +132,15 @@
             thisData = data(:, n);
             violins(n) = Violin(thisData, n, varargin{:});
         end
-        set(gca, 'XTick', 1:size(data, 2));
+        set(AxisHandle, 'XTick', 1:size(data, 2));
         if hascategories && length(cats) == size(data, 2)
-            set(gca, 'XTickLabels', cats);
+            set(AxisHandle, 'XTickLabel', cats);
         end
 
     end
+    
+    if nargout > 0
+        varargout{1} = violins;
+    end
 
 end