diff --git a/.github/workflows/run_book_tests.yml b/.github/workflows/run_book_tests.yml
index 0880ab5..6a2f660 100644
--- a/.github/workflows/run_book_tests.yml
+++ b/.github/workflows/run_book_tests.yml
@@ -15,9 +15,27 @@ jobs:
       - name: Check out RVC3-MATLAB repository (for unit tests)
         uses: actions/checkout@v3
       - run: echo "The ${{ github.repository }} repository has been cloned to the runner."          
+
+      - name: Start display server
+        # Start a display server so Java Swing is available for testing
+        run: |
+          sudo apt-get install xvfb
+          Xvfb :99 &
+          echo "DISPLAY=:99" >> $GITHUB_ENV
+          
       - name: Set up MATLAB
         # This uses the latest release of MATLAB. Can specify "release" if needed.
-        uses: matlab-actions/setup-matlab@v1
+        # Use v2-beta for Java Swing access
+        uses: matlab-actions/setup-matlab@v2-beta
+        with:
+          release: R2023b
+          products: MATLAB Simulink Robotics_System_Toolbox Automated_Driving_Toolbox Computer_Vision_Toolbox Image_Processing_Toolbox Deep_Learning_Toolbox Model_Predictive_Control_Toolbox Optimization_Toolbox ROS_Toolbox Signal_Processing_Toolbox Statistics_and_Machine_Learning_Toolbox Symbolic_Math_Toolbox UAV_Toolbox Control_System_Toolbox
+
+      - name: Print out ver details
+        uses: matlab-actions/run-command@v1
+        with:
+          command: ver
+        
       - name: Run MATLAB Tests
         # Only run tests in folder test/book
         uses: matlab-actions/run-tests@v1
diff --git a/.github/workflows/run_toolbox_tests.yml b/.github/workflows/run_toolbox_tests.yml
index 7b276fe..ab551a5 100644
--- a/.github/workflows/run_toolbox_tests.yml
+++ b/.github/workflows/run_toolbox_tests.yml
@@ -17,6 +17,12 @@ jobs:
       - name: Set up MATLAB
         # This uses the latest release of MATLAB. Can specify "release" if needed.
         uses: matlab-actions/setup-matlab@v1
+      
+      - name: Print out ver details
+        uses: matlab-actions/run-command@v1
+        with:
+          command: ver
+          
       - name: Run MATLAB Tests
         # Only run tests in test/toolbox folder
         uses: matlab-actions/run-tests@v1
diff --git a/test/toolbox/RunAllTests.m b/test/toolbox/RunAllTests.m
new file mode 100644
index 0000000..b834ec1
--- /dev/null
+++ b/test/toolbox/RunAllTests.m
@@ -0,0 +1,28 @@
+% run the unit tests locally
+%
+% run from this folder, will produce test report + coverage.xml file
+
+%% set up the test runner
+import matlab.unittest.plugins.CodeCoveragePlugin
+import matlab.unittest.plugins.codecoverage.CoberturaFormat
+import matlab.unittest.TestRunner
+import matlab.unittest.plugins.codecoverage.CoverageReport
+
+suite = testsuite(".",IncludeSubfolders=false);
+runner = TestRunner.withTextOutput;
+
+% add a coverage report
+reportFile = 'coverage.xml';
+% reportFormat = CoberturaFormat(reportFile);  % XML report
+reportFormat = CoverageReport(".");            % HTML report
+
+plugin = CodeCoveragePlugin.forFolder("../../toolbox", IncludingSubfolders=true, Producing=reportFormat);
+runner.addPlugin(plugin);
+
+
+%% Run all unit tests in this folder
+fprintf('---------------------------------- Run the unit tests ------------------------------------\n')
+
+results = runner.run(suite);
+open(fullfile(".","index.html"))
+
diff --git a/toolbox/@CentralCamera/visjac_p_polar.m b/toolbox/@CentralCamera/visjac_p_polar.m
index c4bf8b0..586b4f7 100644
--- a/toolbox/@CentralCamera/visjac_p_polar.m
+++ b/toolbox/@CentralCamera/visjac_p_polar.m
@@ -20,14 +20,14 @@
 
 function J = visjac_p_polar(cam, rt, Z)
 
-    if numcols(rt) > 1
+    if numel(rt) > 2
         J = [];
         if length(Z) == 1
             % if depth is a scalar, assume same for all points
-            Z = repmat(Z, 1, numcols(rt));
+            Z = repmat(Z, 1, size(rt, 2));
         end
         % recurse for each point
-        for i=1:numcols(rt)
+        for i=1:size(rt, 2)
             J = [J; visjac_p_polar(cam, rt(:,i), Z(i))];
         end
         return;
diff --git a/toolbox/NOTNEEDED/printline.m b/toolbox/NOTNEEDED/printline.m
index cae12a5..d0e7894 100644
--- a/toolbox/NOTNEEDED/printline.m
+++ b/toolbox/NOTNEEDED/printline.m
@@ -40,8 +40,16 @@
 
 % Copyright 2022-2023 Peter Corke, Witold Jachimczyk, Remo Pillat
 
-function out = printline(X, varargin)
-
+function out = printline(X, options)
+    arguments
+        X
+        options.mode (1,1) string = "rpy";
+        options.fmt (1,1) string = ""
+        options.unit (1,1) string {mustBeMember(options.unit, ["rad", "deg"])} = "rad"
+        options.fid (1,1) {mustBeInteger} = 1
+        options.label (1,1) string = ""
+    end
+    
     opt.fmt = [];
     opt.mode = {'rpy', 'euler', 'axang'};
     opt.unit = 'rad';
diff --git a/toolbox/ccxyz.m b/toolbox/ccxyz.m
index 9c2ee92..a8a6f71 100644
--- a/toolbox/ccxyz.m
+++ b/toolbox/ccxyz.m
@@ -19,9 +19,9 @@
 
 function [x,y] = ccxyz(lambda, e)
     xyz = cmfxyz(lambda);
-    if nargin == 1,
+    if nargin == 1
         cc = xyz ./ (sum(xyz')'*ones(1,3));
-    elseif nargin == 2,
+    elseif nargin == 2
         xyz = xyz .* (e(:)*ones(1,3));
         xyz = sum(xyz);
         cc = xyz ./ (sum(xyz')'*ones(1,3));
@@ -33,3 +33,4 @@
         x = cc(:,1);
         y = cc(:,2);
     end
+end
\ No newline at end of file
diff --git a/toolbox/plotellipse.m b/toolbox/plotellipse.m
index c557ccc..8b15166 100644
--- a/toolbox/plotellipse.m
+++ b/toolbox/plotellipse.m
@@ -1,4 +1,3 @@
-function handles = plotellipse(E, varargin)
 %PLOTELLIPSE Draw an ellipse
 %
 %   PLOTELLIPSE(E) draws an ellipse defined by X'EX = 1
@@ -16,8 +15,12 @@
 %   
 %   PLOTELLIPSE(..., Name=Value) specifies additional
 %   options using one or more name-value pair arguments.
-%   Specify the options after all other input arguments.
-% 
+%
+% The ellipse is defined by x' * E * x = s^2 where x is in R^2
+% and s is the scale factor.
+%
+% Options:
+%
 %   confidence - Confidence interval, range 0 to 1. If a confidence interval 
 %                is given then E is interpretted as an inverse covariance
 %                matrix and the ellipse size is computed using an inverse 
@@ -38,8 +41,6 @@
 %                Default: 1
 %
 %
-% The ellipse is defined by x' * E * x = s^2 where x is in R^2
-% and s is the scale factor.
 % For some common cases we require inv(E), for example
 %     - for robot manipulability
 %       \nu inv(J*J') \nu
@@ -67,52 +68,44 @@
 %      % Draw 95% confidence ellipse
 %      PLOTELLIPSE(COVAR, confidence=0.95);
 %
+% References:
+% - Robotics, Vision & Control: Fundamental algorithms in MATLAB, 3rd Ed.
+%   P.Corke, W.Jachimczyk, R.Pillat, Springer 2023.
+%   Appendix C.1.4
 %
 % See also PLOTCIRCLE, PLOT_BOX, CHI2INV.
 
 % Copyright 2022-2023 Peter Corke, Witek Jachimczyk, Remo Pillat
 
-    opt.fillcolor = 'none';
-    opt.alpha = 1;
-    opt.edgecolor = 'k';
-    opt.alter = [];
-    opt.npoints = 40;
-    opt.confidence = [];
-    opt.inverted = 0;
-
-    [opt,arglist,ls] = tb_optparse(opt, varargin);
-
-    % process some arguments
-
-    if ~isempty(ls)
-        opt.edgecolor = ls{1};
-    end
+function handles = plotellipse(varargin)
+    ip = inputParser();
+    ip.KeepUnmatched = true;
+    ip.addRequired("E", @(x) isnumeric(x) && isreal(x));
+    ip.addOptional("center", [0 0], @(x) isnumeric(x) && isreal(x) && (numel(x) == 2));
+    ip.addOptional("ls", "", @(x) ischar(x) || isstring(x));
+    ip.addParameter("fillcolor", "none");
+    ip.addParameter("alpha", 1, @(x) isnumeric(x) && isreal(x) && isscalar(x));
+    ip.addParameter("edgecolor", "k");
+    ip.addParameter("alter", [], @(x) ishandle(x));
+    ip.addParameter("npoints", 40);
+    ip.addParameter("confidence", [], @(x) isnumeric(x));
+    ip.addParameter("inverted", false, @(x) islogical(x));
+    ip.parse(varargin{:});
+    args = ip.Results;
+    arglist = namedargs2cell(ip.Unmatched);
+    E = args.E;
 
     % process the probability
-    if isempty(opt.confidence)
+    if isempty(args.confidence)
         s = 1;
     else
         if exist('chi2inv', 'file') == 2
-            s = sqrt(chi2inv(opt.confidence, 2));
+            s = sqrt(chi2inv(args.confidence, 2));
         else
-            s = sqrt(chi2inv_rvc(opt.confidence, 2));
+            s = sqrt(chi2inv_rvc(args.confidence, 2));
         end
     end
 
-    if ~isempty(arglist) && isnumeric(arglist{1})
-        % ellipse center is provided
-        center = arglist{1};
-        arglist = arglist(2:end);
-    else
-        % default to origin
-        center = zeros(1, size(E,1));
-    end
-
-    % check the ellipse to be altered
-    if ~isempty(opt.alter) && ~ishandle(opt.alter)
-        error('SMTB:plotellipse:badarg', 'argument to alter must be a valid graphic object handle');
-    end
-
     holdon = ishold();
     hold on
 
@@ -125,70 +118,68 @@
         E = [A C;C B];
 
     elseif all(size(E) == [2 2])
-        if ~opt.inverted
+        if ~args.inverted
             E = inv(E);
         end
     else
-        error('ellipse is defined by a 2x2  matrix');
+        error('RVC3:plotellipse:badarg', ...
+            'ellipse is defined by a 1x3 or 2x2  matrix');
     end
 
 
     %% plot an ellipse
 
     % define points on a unit circle
-    th = linspace(0, 2*pi, opt.npoints);
+    th = linspace(0, 2*pi, args.npoints);
     pc = [cos(th);sin(th)];
 
     % warp it into the ellipse
     pe = sqrtm(E)*pc * s;
 
     % offset it to optional non-zero center point
-    center = center(:);
-    if nargin > 1
-        pe = bsxfun(@plus, center(1:2), pe);
-    end
+    pe = args.center(:) + pe;
     x = pe(1,:); y = pe(2,:);
 
     % plot 2D data
 
-    if length(center) > 2
+    if length(args.center) > 2
         % plot 3D data
-        z = ones(size(x))*center(3);
+        z = ones(size(x))*args.center(3);
     else
         z = zeros(size(x));
     end
 
 
-    if strcmpi(opt.fillcolor, 'none')
+    if strcmpi(args.fillcolor, "none")
         % outline only, draw a line
 
         if isempty(ls)
-            if ~isempty(opt.edgecolor)
-                arglist = ['Color', opt.edgecolor, arglist];
+            if ~isempty(args.edgecolor)
+                arglist = [{"Color"}, args.edgecolor, arglist];
             end
         else
-            arglist = [ls arglist];
+            arglist = [{args.ls}, arglist];
         end
 
-        if isempty(opt.alter)
+        if isempty(args.alter)
             h = plot3(x', y', z', arglist{:});
         else
-            set(opt.alter, 'xdata', x, 'ydata', y);
+            set(args.alter, "xdata", x, "ydata", y);
         end
     else
         % fillcolored, use a patch
 
-        if ~isempty(opt.edgecolor)
-            arglist = ['EdgeColor', opt.edgecolor, arglist];
+        if ~isempty(args.edgecolor)
+            arglist = ["EdgeColor", args.edgecolor, arglist];
         end
 
-        arglist = [ls, 'FaceAlpha', opt.alpha, arglist];
+        arglist = [ls, "FaceAlpha", args.alpha, arglist];
 
 
-        if isempty(opt.alter)
-            h = patch(x', y', z', opt.fillcolor, arglist{:});
+        if isempty(args.alter)
+            h = patch(x', y', z', args.fillcolor, arglist{:});
         else
-            set(opt.alter, 'xdata', x, 'ydata', y);
+            set(args.alter, "xdata", x, "ydata", y);
         end
 
     end
diff --git a/toolbox/plotellipsoid.m b/toolbox/plotellipsoid.m
index ef686e5..b6b003c 100644
--- a/toolbox/plotellipsoid.m
+++ b/toolbox/plotellipsoid.m
@@ -10,8 +10,19 @@
 %   
 %   PLOTELLIPSOID(..., Name=Value) specifies additional
 %   options using one or more name-value pair arguments.
-%   Specify the options after all other input arguments.
-% 
+%
+% The ellipsoid is defined by x' * E * x = 1 where x is in R^3.
+%
+% For some common cases we require inv(E), for example
+%     - for robot manipulability
+%       \nu inv(J*J') \nu
+%     - a covariance matrix
+%       (x - \mu)' inv(P) (x - \mu)
+%  so to avoid inverting E twice to compute the ellipsoid, we flag that
+%  the inverse is provided using "inverted".
+%
+% Options:
+%
 %   inverted   - If true, E is inverted (e.g. covariance matrix)
 %                Default: false   
 %   alter      - Alter existing ellipsoids with handle H. This can be used to
@@ -24,124 +35,111 @@
 %                Default: "none"
 %   alpha      - Transparency of the fillcolored ellipsoid: 0=transparent, 1=solid
 %                Default: 1
-%   shadow     - If true, shows shadows on the 3 walls of the plot box
+%   shadow     - If true, shows shadows on the 3 walls of the plot box,
+%                eg. as shown in Fig 8.8 (p 341) of RVC3.
 %                Default: false
 %
+% Other options can be passed to the underlying mesh function, for example:
 %
-% The ellipsoid is defined by x' * E * x = s^2 where x is in R^2
-% and s is the scale factor.
-% For some common cases we require inv(E), for example
-%     - for robot manipulability
-%       \nu inv(J*J') \nu
-%     - a covariance matrix
-%       (x - \mu)' inv(P) (x - \mu)
-%  so to avoid inverting E twice to compute the ellipsoid, we flag that
-%  the inverse is provided using "inverted".
+%     plotellipsoid(diag([1 2 3]), LineStyle=":")
 %
 % - For an unfilled ellipse:
-%   - any standard MATLAB LineStyle such as 'r' or 'b---'.
+%   - any standard MATLAB LineStyle such as "r" or "b---".
 %   - any MATLAB LineProperty options can be given such as 'LineWidth', 2.
 % - For a filled ellipse any MATLAB PatchProperty options can be given
 %
-% Notes::
-% - The 'alter' option can be used to create a smooth animation.
+% Notes:
+% - The alter option can be used to create a smooth animation.
 % - The ellipse is added to the current plot irrespective of hold status.
 %
-% See also PLOTELLIPSE, PLOTCIRCLE, PLOT_BOX, PLOTPOLY, CHI2INV.
+% References:
+% - Robotics, Vision & Control: Fundamental algorithms in MATLAB, 3rd Ed.
+%   P.Corke, W.Jachimczyk, R.Pillat, Springer 2023.
+%   Appendix C.1.4
+%
+% See also PLOTELLIPSE, PLOTSPHERE.
 
 % Copyright 2022-2023 Peter Corke, Witold Jachimczyk, Remo Pillat
 
-function handles = plotellipsoid(E, varargin)
-
-assert(size(E,1) == size(E,2), 'ellipse is defined by a square matrix');
-assert( size(E,1) == 3, 'E must be 3x3 for an ellipsoid');
-
-opt.fillcolor = 'none';
-opt.alpha = 1;
-opt.edgecolor = 'k';
-opt.alter = [];
-opt.npoints = 40;
-opt.shadow = 0;
-opt.inverted = 0;
-
-[opt,arglist,ls] = tb_optparse(opt, varargin);
-
-% process some arguments
-
-if ~isempty(ls)
-    opt.edgecolor = ls{1};
-end
-
-
-if ~isempty(arglist) && isnumeric(arglist{1})
-    % ellipse center is provided
-    center = arglist{1};
-    arglist = arglist(2:end);
-else
-    % default to origin
-    center = zeros(1, size(E,1));
-end
-
-% check the ellipse to be altered
-if ~isempty(opt.alter) && ~ishandle(opt.alter)
-    error('SMTB:plotellipse:badarg', 'argument to alter must be a valid graphic object handle');
-end
-
-if ~opt.inverted
-    E = inv(E);
-end
-
-holdon = ishold();
-hold on
-
-%% plot an ellipsoid
-
-% define mesh points on the surface of a unit sphere
-[Xs,Ys,Zs] = sphere();
-ps = [Xs(:) Ys(:) Zs(:)]';
-
-% warp it into the ellipsoid
-pe = sqrtm(E) * ps;
-
-% offset it to optional non-zero center point
-if nargin > 1
-    pe = bsxfun(@plus, center(:), pe);
-end
-
-% put back to mesh format
-Xe = reshape(pe(1,:), size(Xs));
-Ye = reshape(pe(2,:), size(Ys));
-Ze = reshape(pe(3,:), size(Zs));
-
-
-if isempty(opt.alter)
-    % plot it
-    %             Ce = ones(size(Xe));
-    %             Ce = cat(3, Ce*0.8, Ce*0.4, Ce*0.4);
-    h = mesh(Xe, Ye, Ze, 'FaceColor', opt.fillcolor, ...
-        'FaceAlpha', opt.alpha, 'EdgeColor', opt.edgecolor, arglist{:});
-else
-    % update an existing plot
-    set(opt.alter, 'xdata', Xe, 'ydata', Ye, 'zdata', Ze,  ...
-        arglist{:});
-end
-
-% draw the shadow
-if opt.shadow
-    I = ones(size(Xe));
-    a = [xlim ylim zlim];
-    mesh(a(1)*I, Ye, Ze, 'FaceColor', 0.7*[1 1 1], 'EdgeColor', 'none', 'FaceAlpha', 0.5);
-    mesh(Xe, a(3)*I, Ze, 'FaceColor', 0.7*[1 1 1], 'EdgeColor', 'none', 'FaceAlpha', 0.5);
-    mesh(Xe, Ye, a(5)*I, 'FaceColor', 0.7*[1 1 1], 'EdgeColor', 'none', 'FaceAlpha', 0.5);
-end
-
-
-
-if ~holdon
-    hold off
-end
-
-if nargout > 0
-    handles = h;
-end
+function handles = plotellipsoid(varargin)
+    ip = inputParser();
+    ip.KeepUnmatched = true;
+    ip.addRequired("E", @(x) isnumeric(x) && isreal(x) && all(size(x) == [3 3]));
+    ip.addOptional("center", [0 0 0], @(x) isnumeric(x) && isreal(x) && (numel(x) == 3));
+    ip.addParameter("fillcolor", "none");
+    ip.addParameter("alpha", 1, @(x) isnumeric(x) && isreal(x) && isscalar(x));
+    ip.addParameter("edgecolor", "k");
+    ip.addParameter("alter", [], @(x) ishandle(x));
+    ip.addParameter("npoints", 40);
+    ip.addParameter("shadow", false, @(x) islogical(x));
+    ip.addParameter("inverted", false, @(x) islogical(x));
+    ip.parse(varargin{:});
+    args = ip.Results;
+    arglist = namedargs2cell(ip.Unmatched);
+    E = args.E;
+
+    if ~args.inverted
+        E = inv(E);
+    end
+    
+    holdon = ishold();
+    hold on
+    
+    %% plot an ellipsoid
+    
+    % define mesh points on the surface of a unit sphere
+    [Xs,Ys,Zs] = sphere();
+    ps = [Xs(:) Ys(:) Zs(:)]';
+    
+    % warp it into the ellipsoid
+    pe = sqrtm(E) * ps;
+    
+    % offset it to optional non-zero center point
+    pe = args.center(:) + pe;
+    
+    % put back to mesh format
+    Xe = reshape(pe(1,:), size(Xs));
+    Ye = reshape(pe(2,:), size(Ys));
+    Ze = reshape(pe(3,:), size(Zs));
+    
+    if isempty(args.alter)
+        % plot it
+        %             Ce = ones(size(Xe));
+        %             Ce = cat(3, Ce*0.8, Ce*0.4, Ce*0.4);
+        h = mesh(Xe, Ye, Ze, ...
+            'FaceColor', args.fillcolor, ...
+            'FaceAlpha', args.alpha, ...
+            'EdgeColor', args.edgecolor, ...
+            arglist{:});
+    else
+        % update an existing plot
+        set(opt.alter, 'xdata', Xe, 'ydata', Ye, 'zdata', Ze,  ...
+            arglist{:});
+    end
+    
+    % draw the shadow
+    if args.shadow
+        I = ones(size(Xe));
+        a = [xlim ylim zlim];
+        mesh(a(1)*I, Ye, Ze, ...
+            'FaceColor', 0.7*[1 1 1], ...
+            'EdgeColor', 'none', ...
+            'FaceAlpha', 0.5);
+        mesh(Xe, a(3)*I, Ze, ...
+            'FaceColor', 0.7*[1 1 1], ...
+            'EdgeColor', 'none', ...
+            'FaceAlpha', 0.5);
+        mesh(Xe, Ye, a(5)*I, ...
+            'FaceColor', 0.7*[1 1 1], ...
+            'EdgeColor', 'none', ...
+            'FaceAlpha', 0.5);
+    end
+    
+    if ~holdon
+        hold off
+    end
+    
+    if nargout > 0
+        handles = h;
+    end
 end
diff --git a/toolbox/plotsphere.m b/toolbox/plotsphere.m
index 95dcdee..5615048 100644
--- a/toolbox/plotsphere.m
+++ b/toolbox/plotsphere.m
@@ -1,51 +1,51 @@
 %PLOTSPHERE Draw sphere
 %
-% PLOTSPHERE(C, R, LS) draws spheres in the current plot.  C is the 
-% center of the sphere (1x3), R is the radius and LS is an optional MATLAB 
-% ColorSpec, either a letter or a 3-vector.  
+% PLOTSPHERE(C, R) draws spheres in the current plot.  C is the 
+% center of the sphere (1x3), and R is the radius. 
 %
-% PLOTSPHERE(C, R, COLOR, ALPHA) as above but ALPHA specifies the opacity
-% of the sphere where 0 is transparant and 1 is opaque.  The default is 1.
+% PLOTSPHERE(C, R, COLOR) as above and COLOR is an optional MATLAB 
+% ColorSpec, either a letter or a 3-vector.  
 %
-% If C (Nx3) then N sphhere are drawn and H is Nx1.  If R (1x1) then all
+% If C (Nx3) then N spheres are drawn.  If R (1x1) then all
 % spheres have the same radius or else R (1xN) to specify the radius of
 % each sphere.
 %
 % H = PLOTSPHERE(...) as above but returns the handle(s) for the
 % spheres.
 %
-% Notes::
+% Notes:
 % - The sphere is always added, irrespective of figure hold state.
-% - The number of vertices to draw the sphere is hardwired.
 %
-% Example::
-%         plotsphere( mkgrid(2, 1), .2, 'b'); % Create four spheres
-%         lighting gouraud  % full lighting model
-%         light
+% Example:
+%    plotsphere(mkgrid(2, 1), .2, color="b"); % Create four spheres
+%    lighting gouraud  % full lighting model
+%    light
+%
+% Options:
 %
-% See also: plotpoint, plot_box, plot_circle, plotellipse, plotpoly.
+% color     - surface color as a MATLAB ColorSpec (defaults to "b")
+% alpha     - opaqueness (defaults to 1)
+% edgecolor - edge color as a MATLAB ColorSpec (defaults to "none")
+% n         - number of points on sphere (defaults to 40)
+%
+% See also PLOTELLIPSOID.
 
 % Copyright 2022-2023 Peter Corke, Witold Jachimczyk, Remo Pillat
 
-
-% TODO
-% inconsistant call format compared to other plot_xxx functions.
-
-function out = plotsphere(c, r, varargin)
-
-    opt.color = 'b';
-    opt.alpha = 1;
-    opt.mesh = 'none';
-    opt.n = 40;
-
-    [opt,args] = tb_optparse(opt, varargin);
+function out = plotsphere(c, r, color, opt)
+    arguments
+        c (:,3) {mustBeFloat,mustBeReal};
+        r (1,1) {mustBeFloat,mustBeReal,mustBePositive};
+        color string = "b";
+        opt.color = "";
+        opt.alpha = 1;
+        opt.edgecolor = "none";
+        opt.n = 40;
+    end
     
     % backward compatibility with RVC
-    if length(args) > 0
-        opt.color = args{1};
-    end
-    if length(args) > 1
-        opt.alpha = args{2};
+    if color ~= ""
+        opt.color = color;
     end
     
     daspect([1 1 1])
@@ -71,7 +71,10 @@
         z = r(i)*zs + c(i,3);
                 
         % the following displays a nice smooth sphere with glint!
-        h = surf(x,y,z, ones(size(z)), 'FaceColor', opt.color, 'EdgeColor', opt.mesh, 'FaceAlpha', opt.alpha);
+        h = surf(x,y,z, ones(size(z)), ...
+            'FaceColor', opt.color, ...
+            'EdgeColor', opt.edgecolor, ...
+            'FaceAlpha', opt.alpha);
         % camera patches disappear when shading interp is on
         %h = surfl(x,y,z)
     end
@@ -83,3 +86,4 @@
     if nargout > 0
         out = h;
     end
+end
\ No newline at end of file
diff --git a/toolbox/rpy2jac.m b/toolbox/rpy2jac.m
index 0b110e2..2ed3756 100644
--- a/toolbox/rpy2jac.m
+++ b/toolbox/rpy2jac.m
@@ -1,7 +1,7 @@
 %RPY2JAC Jacobian from RPY angle rates to angular velocity
 %
 % J = RPY2JAC(YPR, OPTIONS) is a Jacobian matrix (3x3) that maps ZYX yaw-pitch-roll angle
-% rates to angular velocity at the operating point RPY=[Y,P,R].
+% rates to angular velocity at the operating point YPR=[Y,P,R].
 %
 % J = RPY2JAC(Y, P, R, OPTIONS) as above but the yaw-pitch-roll angles are passed
 % as separate arguments.
@@ -21,27 +21,30 @@
 
 % Copyright 2022-2023 Peter Corke, Witold Jachimczyk, Remo Pillat
 
-function J = rpy2jac(r, varargin)
+function J = rpy2jac(varargin)
 
-opt.order = {'zyx', 'xyz', 'yxz'};
-[opt,args] = tb_optparse(opt, varargin);
+if numel(varargin{1}) == 3
+    % angle passed as a vector
+    ypr = varargin{1};
+    p = ypr(2);
+    y = ypr(1);
 
-
-% unpack the arguments
-if size(r,2) == 3
-    p = r(:,2);
-    y = r(:,1);
-%     roll = r(:,1);
 elseif nargin >= 3
-    p = args{1};
-    y = r;
-    % roll = args{2}
+    p = varargin{2};
+    y = varargin{1};
+
 else
     error('RVC3:rpy2jac:badarg', 'bad arguments')
 end
 
+if isstring(varargin{end}) || ischar(varargin{end})
+    order = lower(string(varargin{end}));
+else
+    order = "zyx";
+end
+
 
-switch opt.order
+switch order
     case 'xyz'
         J = [
             sin(p)          0       1
diff --git a/toolbox/xaxis.m b/toolbox/xaxis.m
index b40cff2..6a2a165 100644
--- a/toolbox/xaxis.m
+++ b/toolbox/xaxis.m
@@ -12,36 +12,24 @@
 
 % Copyright 2022-2023 Peter Corke, Witold Jachimczyk, Remo Pillat
 
-function xaxis(varargin)
+function xaxis(a1, a2)
 
-    opt.all = false;
-    [opt,args] = tb_optparse(opt, varargin);
-    
-    if isempty(args)
-        [x,y] = ginput(2);
-        mn = x(1);
-        mx = x(2);
-    elseif length(args) == 1
-        if length(args{1}) == 1
+    if nargin == 0
+        set(gca, 'XLimMode', 'auto')
+        return
+    elseif nargin == 1
+        if length(a1) == 1
             mn = 0;
-            mx = args{1};
-        elseif length(args{1}) == 2
-            mn = args{1}(1);
-            mx = args{1}(2);
+            mx = a1;
+        elseif length(a1) == 2
+            mn = a1(1);
+            mx = a1(2);
         end
-    elseif length(args) == 2
-        mn = args{1};
-        mx = args{2};
+    elseif nargin == 2
+        mn = a1;
+        mx = a2;
     end
+    
+    set(gca, 'XLimMode', 'manual', 'XLim', [mn mx])
 
-    if opt.all
-        for a=get(gcf,'Children')',
-            if strcmp(get(a, 'Type'), 'axes') == 1,
-                set(a, 'XLimMode', 'manual', 'XLim', [mn mx])
-                set(a, 'YLimMode', 'auto')
-            end
-        end
-    else
-        set(gca, 'XLimMode', 'manual', 'XLim', [mn mx])
-        set(gca, 'YLimMode', 'auto')
-    end
\ No newline at end of file
+end
\ No newline at end of file
diff --git a/toolbox/xplot.m b/toolbox/xplot.m
index d52d685..84fa3a0 100644
--- a/toolbox/xplot.m
+++ b/toolbox/xplot.m
@@ -9,32 +9,42 @@
 % XPLOT(T, Q) as above but displays the joint angle trajectory versus time
 % given the time vector T (Mx1).
 %
+% Options:
+%
+%  unwrap       - if true, values are considered angles and unwrapped
+%
 % See also JTRAJ, PLOTP, PLOT.
 
 % Copyright 2022-2023 Peter Corke, Witold Jachimczyk, Remo Pillat
 
-function xplot(t, q, varargin)
-if nargin < 2
-    q = t;
-    t = (1:size(q,1))';
-end
-
-opt.unwrap = false;
-opt = tb_optparse(opt, varargin);
-
-if opt.unwrap
-    q = unwrap(q);
-end
+function xplot(varargin)
 
-%clf
-hold on
-plot(t, q(:,1:3))
-plot(t, q(:,4:6), '--')
-grid on
-xlabel('Time (s)')
-ylabel('Joint coordinates (rad,m)')
-legend('q1', 'q2', 'q3', 'q4', 'q5', 'q6');
-hold off
+    ip = inputParser();
+    ip.addOptional("t", []);
+    ip.addRequired("q");
+    ip.addParameter("unwrap", false, @(x) islogical(x));
+    ip.parse(varargin{:});
+    args = ip.Results;
 
-xlim([t(1), t(end)]);
+    q = args.q
+    t = args.t;
+    if isempty(t)
+        t = [0:size(q, 1)-1]';
+    end
+    
+    if args.unwrap
+        q = unwrap(q);
+    end
+    
+    %clf
+    hold on
+    plot(t, q(:,1:3))
+    plot(t, q(:,4:6), '--')
+    grid on
+    xlabel('Time (s)')
+    ylabel('Joint coordinates (rad,m)')
+    legend('q1', 'q2', 'q3', 'q4', 'q5', 'q6');
+    hold off
+    
+    xlim([t(1), t(end)]);
 end
diff --git a/toolbox/yaxis.m b/toolbox/yaxis.m
index d227fba..417717d 100644
--- a/toolbox/yaxis.m
+++ b/toolbox/yaxis.m
@@ -13,20 +13,21 @@
 % Copyright 2022-2023 Peter Corke, Witold Jachimczyk, Remo Pillat
 
 function yaxis(a1, a2)
-	if nargin == 0,
+	if nargin == 0
 		set(gca, 'YLimMode', 'auto')
 		return
-	elseif nargin == 1,
-		if length(a1) == 1,
+	elseif nargin == 1
+		if length(a1) == 1
 			mn = 0;
 			mx = a1;
-		elseif length(a1) == 2,
+		elseif length(a1) == 2
 			mn = a1(1);
 			mx = a1(2);
 		end
-	elseif nargin == 2,
+	elseif nargin == 2
 		mn = a1;
 		mx = a2;
 	end
 
 	set(gca, 'YLimMode', 'manual', 'YLim', [mn mx])
+end
\ No newline at end of file