xold=0.1
T=10000;
st=300;
liap=[0 0];


for r=3:0.005:4
   
l=0;
for i=1:T
   
   xnew=r*xold*(1-xold);
   xold=xnew;
   
   if i>st
       
   der=log(abs(r-2*r*xnew));
   
   l=l+der/(T-st);
   end
   
end
liap=[liap
            r l];
        
plot(liap(2:end,1),liap(2:end,2))
grid on
ylabel('Liapunov exponent','FontSize', 24)
xlabel('r','FontSize', 24)

end