[t,lc,iso]=isochron
lc =
0.4557 0.5550 0.5043
0.4543 0.5547 0.5055
0.4530 0.5544 0.5066
0.4516 0.5541 0.5078
0.4503 0.5538 0.5089
0.4489 0.5535 0.5100
0.4475 0.5531 0.5111
0.4462 0.5528 0.5122
0.4448 0.5524 0.5133
0.4434 0.5520 0.5144
iso =
0.1239 0.1241 0.1244 0.1247 0.1249 0.1252 0.1255 0.1257 0.1260 0.1263 0.1266 0.1269 0.1272 0.1275 0.1278 0.1281 0.1285 0.1288 0.1291 0.1294 0.1298 0.1301 0.1305 0.1308 0.1312 0.1315 0.1319 0.1323 0.1327 0.1330
0.6945 0.6940 0.6936 0.6931 0.6926 0.6922 0.6917 0.6913 0.6908 0.6903 0.6899 0.6894 0.6890 0.6885 0.6880 0.6876 0.6871 0.6867 0.6862 0.6857 0.6853 0.6848 0.6844 0.6839 0.6835 0.6830 0.6825 0.6821 0.6816 0.6812
-0.3398 -0.3346 -0.3293 -0.3241 -0.3188 -0.3136 -0.3085 -0.3033 -0.2981 -0.2930 -0.2878 -0.2827 -0.2776 -0.2725 -0.2674 -0.2624 -0.2574 -0.2523 -0.2473 -0.2423 -0.2374 -0.2324 -0.2274 -0.2225 -0.2176 -0.2127 -0.2078 -0.2029 -0.1980 -0.1931
function [t,lc,iso]=isochron()
alpha1 = 0.1;alpha2 = 3;alpha3 = 3;beta1 = 3;beta2 = 1;beta3 = 1;K1 = 0.5;K2 = 0.5;K3 = 0.5;n1 = 8;n2 = 8;n3 = 8;
system = @(t,Y)[alpha1 - beta1*Y(1)*(Y(3)^n1/(K1^n1+Y(3)^n1));
alpha2*(1-Y(2))*(Y(1)^n2/(K2^n2+Y(1)^n2))-beta2*Y(2);
alpha3*(1-Y(3))*(Y(2)^n3/(K3^n3+Y(2)^n3))-beta3*Y(3)];
x0=[0.455676;0.555023;0.504329];
[t,lc] = ode45(system,0:tau:T,x0);
center = (max(lc)+min(lc))'/2;
iso=[x0-m^0.5*dx, x0+m^0.5*dx];
iso(:,i)=iso(:,i)-tau*feval(system,t,iso(:,i));
if any(abs(iso(:,i)-center)>1.2*m*dx)
iso = [iso(:,1:i-1),iso(:,i+1:end)];
d=sqrt(sum(((iso(:,i)-iso(:,i+1))./dx).^2));
iso = [iso(:,1:i), (iso(:,i)+iso(:,i+1))/2 ,iso(:,i+1:end)];
iso = [iso(:,1:i), iso(:,1:i)];
if (mod(-t,T)<=tau/2) && (-t<k*T+tau)
plot3(lc(:,1),lc(:,2),lc(:,3),'r');
if min(abs(mod(-t,T)-phases))<tau/2
plot3(iso(1,:),iso(2,:),iso(3,:),'k-');