Solving a system of Non Linear Differential Equations

Hello, I want to solve the following system of non-linear differential equations numerically. Please provide guidance.

[KSSV]

Have a look on ode45.

[kdv0]

Hey I did this. Is it correct? The graphs are straight lines...

x0 = 0;
y0 = 0;
z0 = 10;
Y0 = [x0 y0 z0];
tspan = [0 1000];
[t, Y] = ode45(@rate, tspan, Y0);

Warning: Failure at t=9.804721e-03. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.775558e-17) at time t.

x = Y(:,1);
y = Y(:,2);
z = Y(:,3);
plot(t, x)

function dYdt = rate(~,Y)

        r=10;
        a=2;
        b=2;
        g=1;
        m=0.5;
        c=1;
        p=1;
        q=4;
        s=1;
        x = Y(1);
        y = Y(2);
        z = Y(3);
                
        dxdt = r*x*(1-x/z)-a*m*x*y/(1+g*m*x);
        dydt = b*m*x*y/(1+g*m*x)-c*y;
        dzdt = p*z^2+q*z+s;
        
        dYdt = [dxdt;
                dydt;
                dzdt];
end

Answer 1

[Sam Chak]

Hi @kdv0

Apart from the incorrect initial value for z, which should be , the rest of the information in the code is correct. The response for x behaves as expected since it starts from the equilibrium point.

Both time derivatives for x and y are zero at the beginning because the states x and y start from the equilibrium point.

Even though the time derivative for z is decoupled from the influence of the variations in x and y, the rate of change exhibits a quadratic positive behavior when z is greater than zero. This is attributed to the positive values of the coefficients p, q, and s.

In simpler terms, the rate of change increases rapidly as z moves away from zero in the positive direction. Consequently, the response of z becomes unstable, leading to an explosive behavior and integration failure at time .

p   = 1;
q   = 4;
s   = 1;
z   = linspace(-14, 10, 2001);
dz  = p*z.^2 + q*z + s;

plot(z, dz), grid on, xlabel('z'), ylabel('dz'), title('dz/dt')
xline(0, '-'), yline(0, '-')

Answer 2

[KSSV]

This is the code..I am getting hight values...you may check and modify the code.

tspan = [0 1];
y0 = [0 0 100];
sol = ode45(@odefun,tspan,y0) ; 

Warning: Failure at t=9.804721e-03. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.775558e-17) at time t.

plot(sol.x,sol.y)

function dydt = odefun(t,y)
r = 10; 
a = 2 ;
b = 2 ; 
g = 1 ;
m = 0.5 ;
c = 1 ;
p = 1 ;
s = 1 ;
q = 4 ; 

dydt = zeros(3,1) ; 
dydt(1) = r*y(1)*(1-y(1)/y(3))-a*m*y(1)*y(2)/(1+g*m*y(1)) ;  
dydt(2) = b*m*y(1)*y(2)/(1+g*m*y(1))-c*y(2) ; 
dydt(3) = p*y(3)^2+q*y(3)+s ; 
end