How to plot isochrons in 3d

I have code for plotting isochrons for a system of ODE's in 2 dimensions, but I need to plot isochrons for a system of 3 ODEs in 3 dimensions. I have already made some adjustments to the plot commands and such, but it is not running the way I intended. Any help or advice would be greatly appreciated !

[t,lc,iso]=isochron

   0

   0

   0

   0

   0

   0

t = -4.3825

lc = 2001×3

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 = 3×470

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)];
T=1.395*pi;     % period of the oscillator
%x0 = [0.226623;0.174165;0.243654];
x0=[0.455676;0.555023;0.504329];   % initial condition (must be on the limit cycle)
num_chron=5;       % number of isochrons plotted
phases=0:T/num_chron:T;
tau=T/2000;        % time step of integration
m=200;              % spatial grid
k=0;                % the number of skipped cycles
[t,lc] = ode45(system,0:tau:T,x0);  % call the DE
dx=(max(lc)-min(lc))'/m;         % spatial resolution  
center = (max(lc)+min(lc))'/2;   % center of the limit cycle   
iso=[x0-m^0.5*dx, x0+m^0.5*dx];  % isochron?s initial segment
for t=0:-tau:-(k+1)*T            % backward integration
    
    for i=1:size(iso,2)
        iso(:,i)=iso(:,i)-tau*feval(system,t,iso(:,i)); % move one step
    end
    
    i=1;
    while i<=size(iso,2)            % remove infinite solutions
        if any(abs(iso(:,i)-center)>1.2*m*dx)    % check boundaries
            iso = [iso(:,1:i-1),iso(:,i+1:end)]; % remove
        else i=i+1;
        end
    end
    i=1;
    while i<=size(iso,2)-1
        d=sqrt(sum(((iso(:,i)-iso(:,i+1))./dx).^2)); % normalized distance
        if d > 2                    % add a point in the middle
            iso = [iso(:,1:i), (iso(:,i)+iso(:,i+1))/2 ,iso(:,i+1:end)];
        end
        if d < 1             % remove the point
            iso = [iso(:,1:i), iso(:,1:i)];
        else i=i+1;
        end
    end
    if (mod(-t,T)<=tau/2) && (-t<k*T+tau)
        cla;
        plot3(lc(:,1),lc(:,2),lc(:,3),'r'); % plot the limit cycle
        hold on; 
    end
    if min(abs(mod(-t,T)-phases))<tau/2   % plot the isochrons
        plot3(iso(1,:),iso(2,:),iso(3,:),'k-');
        disp(i<=size(iso,2)-1)
        drawnow;
    end        
end
xlim(0:1)
ylim(0:1)
zlim(0:1)
hold off;
end

Answer 1

[Arun]

Hey Braden,

I understand that you want a 3-D plot of “isochrons” for a 3-ODEs and have encountered difficulties in obtaining the desired outcomes.

Here is a way to obtain 3-D plot:

if (mod(-t,T)<=tau/2) && (-t<k*T+tau)
        cla;
        scatter3(lc(:,1),lc(:,2),lc(:,3),'r'); % plot the limit cycle
        hold on; 
end
if min(abs(mod(-t,T)-phases))<tau/2   % plot the isochrons
        scatter3(iso(1,:),iso(2,:),iso(3,:),'k');
        disp(i<=size(iso,2)-1)
        drawnow;
end  

Output plots for your reference:

Hope this helps.