% set initial value of mu
mu=0;
% find a fixed point using Ni iterations of Newton's method 
Ni=6;
u=-0.5*ones(1,Ni);
for i=2:Ni,
    u(i)=u(i-1)-g(u(i-1),mu)/gu(u(i-1));
end;
% find initial tangent vector
nullv=null([gu(u(Ni)) gmu(mu)]);
% number of Newton iterations
NN=5;
% step size
ds=0.06;
% number of steps taken
tot=100;
% initialise variables. dat stores output, ei stores eigenvalue
J=zeros(2);
v=zeros(2,NN);
dat=zeros(2,tot+1);
ei=zeros(1,tot+1);
% record first data point
dat(:,1)=[u(Ni) mu]';
ei(1)=gu(u(Ni));
for j=1:tot, % outer loop
    % record the last point found
    vold=dat(:,j);
    % linear prediction of solution
    v(:,1)=dat(:,j)+ds*nullv;
    for i=2:NN, % inner loop
        f=getf(v(:,i-1),vold,nullv,ds);
        J(1,:)=[gu(v(1,i-1)) gmu(v(2,i-1))];
        J(2,:)=nullv';
        v(:,i)=v(:,i-1)-J\f;
    end;
    % check for convergence
    disp([j norm(f)]);
    % store results of Newton iteration
    dat(:,j+1)=v(:,NN);
    % store stability information
    ei(j+1)=J(1,1);
    % calculate new tangent vector and normalise
    nullv=(dat(:,j+1)-dat(:,j))/ds;
    nullv=nullv/norm(nullv);
end;
% plot results
subplot(2,1,1);
plot(dat(2,:),dat(1,:))
xlabel('\mu'); ylabel('u');
grid;
subplot(2,1,2);
plot(dat(2,:),ei)
xlabel('\mu'); ylabel('eigenvalue');
grid;