function mweode15isparse ()

  a = 0.5;
  tspan = [0; 100 * a];
  x0 = [0; 0.01];
  x0dot = [0.0; 0.0];
  options = odeset ('Jacobian', @J, "RelTol", 1e-4, "AbsTol", 1e-6);
  [t, x] = ode15i (@f, tspan, x0, x0dot, options);

  plot (x(:, 1), x(:, 2));
  xlabel ("x");
  ylabel ("x'");
  title ("Van der Pol Phase Portrait");

  ## x'' - a * (1 - x^2) * x' + x = 0
  ## Let x1 := x and x2 := x' then
  ## [x1', x2'] = [x2, a * (1 - x1^2) * x2 - x1]
  ## 0 = [x2 - x1', a * (1 - x1^2) * x2 - x1 - x2']

  function dxdt = f(t, x, xdot)
    dxdt = [x(2)-xdot(1), a*(1-x(1)^2)*x(2)-x(1)-xdot(2)];
  end

  function [dfdx, dfdxdot] = J(t, x, xdot)
    dfdx = sparse ([0, 1; -2*a*x(1)*x(2)-1, a*(1-x(1)^2)]);
    dfdxdot = sparse ([-1, 0; 0, -1]);
  end

end

