How do I solve differential equation in while/for loop using ODE45?

Hi, I'm trying to simulate Rossler Oscillators with a 4 node system. . I would like to generalize this function by running this in a while/for loop, so that I can solve for 8-node, 12-node, n-node by just changing this in the function.

The initial conditions and ode function, etc. are in the main code. Below is the actual working code for simulating a 4 node system,

function yprime = rossler(t,y,K)

a = 0.18; b = 0.2; c = 5.7;

f = [
-(y(2)+y(3)); %dy1
(y(1) + a*y(2)); %dy2
b+y(3)*(y(1)-c); %dy3

-(y(5)+y(6)); %dy4
(y(4) + a*y(5)); %dy5
b+y(6)*(y(4)-c); %dy6

-(y(8)+y(9)); %dy7
(y(7) + a*y(8)); %dy8
b+y(9)*(y(7)-c); %dy9

-(y(11)+y(12)); %dy10
(y(10) + a*y(11)); %dy11
b+y(12)*(y(10)-c); %dy12
yprime = f - (K*y);

Below is the code I am trying to generalize by using for/while loop. Correct me if I am wrong. Thank you!

function yprime = cycfun(t,y,K)
  a = 0.18; b = 0.2; c = 5.7;
  n = 11; i = 1;
  while i < n
      i = i + 3;
  f = [
  -(y(i+1)+y(i+2)); %dy1
  (y(i) + a*y(i+1)); %dy2
  b+y(i+2)*(y(i)-c);]; %dy3
  end
  yprime = f - (K*y);

When I try running this code, I get matrix dimensions must agree. I know I'm messing it up somewhere, just trying to figure out! Thanks!

Below is the main code where I call the function:

      clear all; clc;
      tspan = [0,500]; 
      options = odeset('stats', 'on');

      %adjacency matrix, diagonal matrix, laplacian matrix
      A = [0 1 1 1; 1 0 1 1; 1 1 0 0; 1 1 0 0]; 
      D = [3 0 0 0; 0 3 0 0; 0 0 2 0; 0 0 0 2];
      L = D - A;

      %feedback structure, coupling factor, kronecker product
      H = [1 0 0; 0 0 0; 0 0 0];
      C = 0.001;
      K = kron(L,H);
      K = C*K;

      %initial condition and ODE function
      yzero = [0.093;0.3;0.065;0.093;0.3;0.065;0.093;0.3;0.065;0.093;0.3;0.065];
      [T,Y] = ode45(@cycfun, tspan, yzero,options, K);

      %plotting square difference
      %--- NODE 1 & 2 ---%
      %Square difference of node 1 and 2
      q1 = ((Y(:,1) - Y(:,4)).^2) + ((Y(:,2) - Y(:,5)).^2) +  ((Y(:,3) - Y(:,6)).^2);
      w1 = q1;

      %Plotting node 1 and 2
      t1 = linspace(0,500,length(w1));
      figure,plot(t1,w1);
      title('Coupled Rossler Oscillator');
      xlabel('t');
      ylabel('(x11 - x21)^2 + (x12 - x22)^2 + (x13 - x23)^2');

      %--- NODE 3 & 4 ---%
      %Square difference of node 3 and 4
      q2 = ((Y(:,7) - Y(:,10)).^2) + ((Y(:,8) - Y(:,11)).^2) +  ((Y(:,9) - Y(:,12)).^2);
      w2 = q2;

      %Plotting node 1 and 2
      t2 = linspace(0,500,length(w2));
      figure,plot(t2,w2);
      title('Coupled Rossler Oscillator');
      xlabel('t');
      ylabel('(x31 - x41)^2 + (x32 - x42)^2 + (x33 - x43)^2');

[Nicolas Schmit]

Please post the error message

[Sreenath Umagandhi]

ERROR: Index exceeds matrix dimensions.

Error in rossler (line 22) -(y(14)+y(15)); %dy13

Error in odearguments (line 90) f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115) odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in cyccode (line 18) [T,Y] = ode45(@rossler, tspan, yzero,options, K);

[KSSV]

I suspect problem is with your input...you should input y whose size should be 12...else error..show us how did you call the function?

[Sreenath Umagandhi]

@KSSV, thank you. I have updated my question with the code where I call the function.

Answer 1

[Torsten]

function yprime = cycfun(t,y,K)
  a = 0.18; b = 0.2; c = 5.7;
  n = 12; i = 1; f = [];
  while i < n 
    f = [f;-(y(i+1)+y(i+2));(y(i) + a*y(i+1)); b+y(i+2)*(y(i)-c)];
    i = i+3;
  end
  yprime = f - (K*y);

Best wishes

Torsten.

[Sreenath Umagandhi]

Thank you! I tried this, but not sure what you are trying to do here. I have updated the question with whole code.

However, I still did get some error:

Subscript indices must either be real positive integers or logicals.

Error in cycfun (line 5) f = [f;-(y(i+1)+y(i+2));(y(i) + a*y(i+1)); b+y(i+2)*(y(i)-c)];

[Torsten]

Code corrected ; it should work now.

Best wishes

Torsten.

[Sreenath Umagandhi]

Torsten,

Thanks a lot! It worked like how I wanted :)