trying to implement nonlinear ODE equation

i have three nonlinear oscillator equation

1 ) dx/dt = (mu-r^2)x - w*y +K*F

2) dy/dt = (mu-r^2)y + w*x

3) dw/dt = (-K*F)*(y/(sqrt(x^2+y^2)))

here,

F is any input signal

K is coupling strength (K > 0 is any value)

w is oscillator intrinsic frequency

mu is constant (mu > 0 any value)

r is taking as a constant here .. r=1

third equation dw/dt is learning rule for intrinsic frequency (as w learn frequency of input signal F) and i am getting stuck in this part of equation that how could i implement this part and integrate it in matlab code ?

x and y are state variables.

the following is matlab code without third equation and i want to integrate third equation in this code.

*******************************************

function dz = myeqd(t,y,ti,fi)

dz = zeros(2,1);

mu=0.7;

r=1;

K=1;

w=40;

F=interp1(ti,fi,t);

dz(1)= (mu - r^2)*y(1) - w*y(2) +K*F;

dz(2) = (mu - r^2)*y(2) + w*y(1);

********************************************

calling this myeqd function

Tspan= [0.1:0.1:10]; % time vector

t= Tspan;

fi = cos(2*pi*Tspan); % perturbation

ti=Tspan;

[T,Y]=ode45(@(t,y) myeqd(t,y,ti,fi),Tspan,[0;1]);

plot (T,Y)

**********************************************

Answer 1

[Azzi Abdelmalek]

The third equation will be

dz(3) = (-K*F)*y(2)/sqrt(y(1)^2+y(2)^2)

with

dz = zeros(3,1);

[Shritej Chavan]

Hi Azzi sir,

I am trying the same thing but for four nonlinear oscillator equation. First three are the same. The fourth one is for learning amplitude of the input signal given by following learning rule

da/dt = eta*x*F where eta - coupling constant, F - input signal, x -state variable

 function dz = myeqd(t,x,fi,ti, epsilon, mu, gamma, eta) 
%gamma - speed  of convergence
%mu - radius of limit cycle 
%epsilon & eta - coupling constant
F_t = interp1(ti,fi, t);
k = gamma*(mu - (x(1)^2 + x(2)^2))*x(1) - x(3)*x(2) + epsilon*(F_t) ;
y = gamma*(mu - (x(1)^2 + x(2)^2))*x(2) + x(3)*x(1) ;
w = -epsilon*(F_t)*x(2)/(x(1)^2 + x(2)^2) ; 
a = eta*x(1)*(F_t);
dz = [k;y;w;a] ; 
end

In the command window

ti = 0:0.1:2000 ;
t = ti ;
fi = sin(40*t) ;
[t,x] = ode45(@(t,x)myeqd(t,x,fi,ti,0.9, 1,17, 0.9), t , [1 0 45 1] ) ;

*The problem is the amplitude does converge to a proper value but the frequency does and the amplitude keeps increasing after time when the oscillator has learned the correct frequency.

I really need your help with this one. Hope you will answer.