clear
hold on
abeg=input(' Enter the value of a to begin: ');
astep=input(' Enter the step in a: ');
afin=input(' Enter the final a: ');
nrenorm=input(' Enter number of steps between renormalisations: ');
nmax=input(' Enter the number of steps per phase point: ');
nlyap=nmax/nrenorm;
a=abeg:astep:afin;
na=length(a);
lyapsum(1:na)=0;
x(1,1:na)=0.501;
for nr=1:nlyap
  len(1,1:na)=1;
  for i=1:nrenorm
    iter=(nr-1)*nrenorm+i;
    x(iter+1,:)=a.*x(iter,:).*(1-x(iter,:));
    len(i+1,:)=len(i,:).*a.*(1-2*x(iter,:));
  end
  for k=1:na
    lyapsum(k)=lyapsum(k)+log(abs(len(nrenorm+1,k)));
  end
end
lyapsum(:)=lyapsum(:)/nmax;
plot(a,lyapsum,'r')