Hello,

I would like to plot gait diagrams (the lines you can see on the picture belowà from the solutions obtained with a NL oscillator (composed with 8 coupled odes). Here the result that I would like to obtain.

Initial plot:

Desired plot

I would like to obtain 4 lines corresponding to the 4 elliptic trajectories obtained with the NL oscillator. The four lines should be done like this. When the trajectory is above 0, the line should be colored in green. When the trajectory is below 0, the line should be colored in black.

**May you help me to define this kind of graph called gait diagrams from the solution of the NL oscillator ?**

Here you can find my maple code:

**K:=Matrix([<0, -1, 1, -1>,<-1, 0, -1, 1>,<-1, 1, 0,-1>,<1, -1, -1,0>]);**

**for i to 4**

**do **

**r[i]:=sqrt((u[i](t))^2+(v[i](t))^2):**

**omega[i]:=omega[sw]/(1+exp(b*v[i](t)))+omega[st]/(1+exp(-b*v[i](t))):**

**Equ[i]:=diff(u[i](t),t)=Au*(1-r[i]^2)*u[i](t)-omega[i]*v[i](t):**

**Eqv[i]:=diff(v[i](t),t)=Av*(1-r[i]^2)*v[i](t)+omega[i]*u[i](t)+MatrixVectorMultiply(K,<seq(v[i](t),i=1..4)>)[i]:**

**EqSys[i]:=[Equ[i],Eqv[i]]:**

**end do:**

**paramsCycle:=omega[st]=4*2*Pi,omega[sw]=2*Pi,Au=5,Av=50,b=100;**

**params:=paramsCycle;**

**Differential system **

**sys:=map(op,eval([seq(EqSys[i],i=1..4)],[params]));**

**ic:=[u[1](0)=0, v[1](0)=0,u[2](0)=0, v[2](0)=-0.1,u[3](0)=0, v[3](0)=0.1,u[4](0)=0, v[4](0)=0.1];**

**Résolution1**

**res:=dsolve([sys[],ic[]],numeric):**

**Initial boundaries**

**tcalc:=4;**

**ic2:=[seq(u[i](0)=eval(u[i](t), res(tcalc)),i=1..4),seq(v[i](0)=eval(v[i](t), res(tcalc)),i=1..4)];**

**Résolution2**

**res:=dsolve([sys[],ic2[]],numeric):**

**tmax:= 40:**

**numpts:=100*tmax:**

**plots:-odeplot(res,[t,v[1](t)],0..tmax,thickness=2, view=[0..5, -1.5..1.5],numpoints = numpts);**

**plots:-odeplot(res,[t,v[2](t)],0..tmax,thickness=2, view=[0..5, -1.5..1.5],numpoints = numpts);**

**plots:-odeplot(res,[t,v[3](t)],0..tmax,thickness=2, view=[0..5, -1.5..1.5],numpoints = numpts);**

**plots:-odeplot(res,[t,v[4](t)],0..tmax,thickness=2, view=[0..5, -1.5..1.5],numpoints = numpts);**

**plots:-odeplot(res,[seq([t,v[i](t)+i*5],i=1..4)],0..tmax,thickness=2,view=[0..5,0..25], numpoints = numpts);**

Thanks a lot for your help