function plot_vector_field(X, Y, W, beta)
% Taken from AIH's mctj_demo.m. Given its own m-file to be available to
% other functions

switch size(X,2)
    case 2
        xlim = get(gca, 'XLim');
        ylim = get(gca, 'YLim');
        xi = linspace(xlim(1), xlim(2), 10);
        yi = linspace(ylim(1), ylim(2), 10);
        [Fx, Fy] = meshgrid(xi, yi);
        F = [Fx(:), Fy(:)];
        F2 = mctlibj.MCTReg.applyVectorField(F, Y, X, W, beta);
        quiver(F(:,1), F(:,2), F(:,1) - F2(:,1), F(:,2) - F2(:,2), 'k');

    case 3
        xlim = get(gca, 'XLim');
        ylim = get(gca, 'YLim');
        zlim = get(gca, 'ZLim');
        %
        xi = linspace(xlim(1), xlim(2), 5);
        yi = linspace(ylim(1), ylim(2), 5);
        zi = linspace(zlim(1), zlim(2), 5);

        [Fx, Fy, Fz] = meshgrid(xi, yi, zi);
        F = [Fx(:), Fy(:), Fz(:)];
        F2 = mctlibj.MCTReg.applyVectorField(F, Y, X, W, beta);
        quiver3(F(:,1), F(:,2), F(:,3), F(:,1) - F2(:,1), F(:,2) - F2(:,2), F(:,3) - F2(:,3), 'k');
end
