Unable to solve stiff differential equation

I have solved many nonlinear differential using ode45 previously. 
For the current stiff nonlinear differential equation, I used ode45, ode15s, ode23s, ode23t and ode23tb but unable to get final results.
I got constant warnings which is mentioned below.

Warning: Failure at t=8.048300e+00.  Unable to meet integration tolerances without reducing the step size below the
smallest value allowed (2.859330e-14) at time t.

Warning: RelTol has been increased to 2.22045e-14. 

The code is below.

clear all;clc;
op=odeset('relTol',1e-14,'absTol',1e-14);
global psi D B H M epsilon N
epsilon=0.01;
M=epsilon*1; 
H=epsilon*1; 
psi=5;    
D=epsilon*1; 
B=epsilon*0.5;    
N=epsilon*0.17;   
w= 0.1:0.01:5;
for i=1:1:length(w)
THETA=w(i);
[t,y]=ode23s(@(t,x)F(t,x,THETA),[0 20],[0 0],op);
end

function dy = F(t,y,THETA) 
global psi D B H M N
dy(1) = y(2);
dy(2) = psi.*cos(THETA.*t)-(1-M.*sin(THETA.*t)).*y(1)-B.*y(2)-(H.*(cos(THETA.*t))).*((y(1)).^2)-(D.*(sin(THETA.*t))-N)*((y(1))^3);
dy=dy';
end

Answer 1

[Sam Chak]

Hi @AJMIT KUMAR,

It appears that the system is inherently unstable for . If you simulate the system for a longer time, then the system will explode. Small value of Θ takes longer to explode.

% op = odeset('relTol', 1e-14, 'absTol', 1e-14);

global psi D B H M epsilon N
epsilon = 0.01;
M   = epsilon*1; 
H   = epsilon*1; 
psi = 5;    
D   = epsilon*1; 
B   = epsilon*0.5;    
N   = epsilon*0.17;   
w   = 0.19:0.001:0.195;       % theta 0.198 will stop integrating at around 20 sec
% w   = zeros(1, 2);          % stable
% w   = 0.01:0.001:0.014;     % stable

for i = 1:1:length(w)
    THETA  = w(i);
    [t, y] = ode45(@(t,x) F(t, x, THETA), [0 20], [0 0]);
    plot(t, y), hold on
end
hold off, grid on, xlabel('t, (sec)')

function dy = F(t, y, THETA)
    global psi D B H M N
    dy(1) = y(2);
    dy(2) = psi*cos(THETA*t) - (1 - M*sin(THETA*t))*y(1) - B*y(2) - H*cos(THETA*t)*y(1)^2 - (D*sin(THETA*t) - N)*y(1)^3;
    dy    = dy';
end

[Sam Chak]

Hi @AJMIT KUMAR

Just to add something; your nonlinear system is unstable

because the coefficient of the cubic growth term is "POSITIVE". In other words, if y gets large, it rapidly "adds" energy to . On the contrary, if you make , then the system will be stable because it rapidly "dissipates" energy in .

If the following example, the sign of N is switched, and stays within .

Note: I didn't perform rigorous mathematical proof to show that the nonlinear system is stable for . It is just my intuitionistic logic.

global psi D B H M epsilon N
epsilon = 0.01;
M   = epsilon*1; 
H   = epsilon*1; 
psi = 5;    
D   = epsilon*1; 
B   = epsilon*0.5;    
N   = - epsilon*0.17;       % <--- this coefficient must be negative    
w   = 0.19:0.001:0.195;     % theta 0.198 will stop integrating at around 20 sec
% w   = zeros(1, 2);        % stable
% w   = 0.01:0.001:0.014;   % stable

for i = 1:1:length(w)
    THETA  = w(i);
    [t, y] = ode45(@(t,x) F(t, x, THETA), [0 3000], [0 0]);
    plot(t, y(:,1)), hold on
end
hold off, grid on, xlabel('t, (sec)'), ylabel('y(t)')
title({'Plot of the solution for $y_{1}(t)$'}, 'interpreter', 'latex', 'fontsize', 16)

idx  = find(t>2000);
ymax = y(:,1);
max(ymax(idx))

ans = 5.2387

function dy = F(t, y, THETA)
    global psi D B H M N
    dy(1) = y(2);
    dy(2) = psi*cos(THETA*t) - (1 - M*sin(THETA*t))*y(1) - B*y(2) - H*cos(THETA*t)*y(1)^2 - (D*sin(THETA*t) - N)*y(1)^3;
    dy    = dy';
end

Answer 2

[Oguz Kaan Hancioglu]

I think matlab solver couldn't solve the nonlinear probrem within your tolerances. You can use fixed step size which is very small dt values (e-5, or e-6) to solve your problem.

https://www.mathworks.com/matlabcentral/answers/92961-how-do-i-use-a-fixed-step-size-with-ode23-and-ode45-in-matlab