function cont01

mu0   = -1;
[t z] = ode45(@(t,x) func(x,mu0), [0 5], 0);                  %#ok<ASGLU>
x0    = z(end);

[xs mus] = run_cont(@func, @func_dx, x0, mu0, 0.1, 40);
plot(mus, xs)

end

function y = func(x,mu)
y = (mu+x-x^3)*( (x-2)^2 - mu+2 );
end

function J = func_dx(x,mu)
a  = mu+x-x^3;
b  = (x-2)^2 - mu+2;
ap = 1-3*x^2;
bp = 2*(x-2);

J = ap*b + a*bp;
end

function [xs mus] = run_cont(f, fx, x0, mu0, h, Nmx)
xs  = x0;
mus = mu0;

  for It=2:Nmx
    [x0 mu0] = step(f,fx,x0,mu0,h);
    xs  = [xs x0];                           %#ok<AGROW>
    mus = [mus mu0];                         %#ok<AGROW>
  end
end

function [x1 mu1] = step(f, fx, x0, mu0, h)
mu1 = mu0+h;
x1  = correct(f, fx, x0, mu1);
end

function x1 = correct(f, fx, x0, mu0)

J  = fx(x0,mu0);
x1 = x0 - J\f(x0,mu0);

while norm(x1-x0)>1.0e-4
  x0 = x1;
  x1 = x1 - J\f(x1, mu0);
end

end
