function cont02

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

[zs ts] = run_cont(@func, @func_dx, @func_dmu, x0, mu0, 0.1, 80); %#ok<NASGU>
plot(zs(2,:), zs(1,:))

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 J = func_dmu(x,mu)
a  = mu+x-x^3;
b  = (x-2)^2 - mu+2;
ap = 1;
bp = -1;

J = ap*b+a*bp;
end

function [zs ts] = run_cont(f, fx, fmu, x0, mu0, h, Nmx)
z0 = [ x0; mu0 ];
zs = z0;
t0 = [0;1];
ts = t0;

  for It=2:Nmx
    [z0 t0] = step(f,fx,fmu,z0,t0,h);
    zs  = [zs z0];                           %#ok<AGROW>
    ts  = [ts t0];                           %#ok<AGROW>
  end
end

function [x1 t1] = step(f, fx, fmu, z0, t0, h)
z1      = z0 + h*t0;
[x1 t1] = correct(@xf, @xfx, z1);

  function y = xf(x)
    y(1,1) = f( x(1), x(2) );
    y(2,1) = t0'*(x-z0) - h;
  end

  function J = xfx(x)
    J(1,1) = fx( x(1), x(2) );
    J(1,2) = fmu( x(1), x(2) );
    J(2,:) = t0';
  end
end

function [x1 t1] = correct(f, fx, x0)

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

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

J  = fx(x1);
t1 = J\[0;1];
t1 = t1/norm(t1);

end
