Atractor de Rossler Matlab

%Resolución del atractor de Rossler mediante un esquema
%Predictor-Corrector de dos pasos, con paso temporal variable. Se utiliza
%el estimador de Milne para controlar el paso temporal dt1.Para ver los
%resultados sin la influencia de las condiciones iniciales plot a partir de
%n=20000.

clear

n=100000;
dt1=0.01;
dt2=0.01;
tol=0.01;
a=0.2;
b=0.2;
c=5.7;
u=zeros(n,3);%AM2
up=zeros(n,3); %AB2
tpc=zeros(n,3); %error de truncación
dt=zeros(n,1); %paso temporal en función del tiempo
t=zeros(n,1); %tiempo acumulado t=n*dt
modtpc=zeros(n,1);%módulo de tpc en función del tiempo
b0c=0.5;
b1c=0.5;
b2p=-dt1/(2*dt2);
b1p=1-b2p;
A=[0 -1 1; 2 a 0; u(1,3) 0 -c];
B=[0 0 b];
I=eye(3);
u(1,:)=[1 1 1];
t(1)=0;
dt(1)=dt1;

%Arranque Euleru(2,:)=u(1,:)*(I+dt1.*A);
t(2)=t(1)+dt1;
dt(2)=dt1;

for i=2:n
A=[0 1 u(i,3); -1 a 0; -1 0 -c];
for j=1:100
b2p=-dt1/(2*dt2);
b1p=1-b2p;up(i+1,:)=u(i,:)+dt1*(b1p*(u(i,:)*A+B)+b2p*(u(i-1,:)*A+B));
c3=-1/12;
c3s=1/6+dt2/(4*dt1);
tpc(i+1,:)=c3*(u(i,:)-up(i+1,:))/(c3s-c3);
modtpc(i+1)=norm(tpc(i+1,:));if modtpc(i+1) > tol
dt1=dt1/2;
else
dt1=dt1*(tol/modtpc(i+1))^(1/3);
end
end
dt2=dt1;
t(i+1)=t(i)+dt1;
dt(i+1)=dt1;u(i+1,:)=u(i,:)+dt1*(b0c*(up(i+1,:)*A+B)+b1c*(u(i,:)*A+B));
end

%plot3(u(:,1),u(:,2),u(:,3))
%plot(t,dt)
%plot(t(1:1000),modtpc(1:1000))

%Estabilidad

p=0.1;[x,y]=meshgrid(-3:p:3,3:-p:-3);
[v]=[1;0.9;0.8;0.7;0.6;0.5;0.4;0.3;0.2];

%Raiz esquema Runge-Kutta
[rhoRK]=abs(1+complex(x,y)+complex(x,y).^2/2);
%contour(x,y,rho,v)

%Raices AB2[r1]=1/2+complex(x,y)*3/4+1/2((1+3/2*complex(x,y)).^2-2*complex(x,y).^(1/2);
[r2]=1/2+complex(x,y)*3/4-1/2((1+3/2*complex(x,y)).^2-2*complex(x,y).^(1/2);
[rhoAB2]=max(abs(r1),abs(r2));

%Raices Predictor-Corrector...