Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN

Question

The following code uses rk4 to simulate the dynamics defined via fe function. However, I encountered the warning Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. .

I felt that the function defined in fe, sometimes the Numerator equals zero, thus creates singularity. But I don't know how to avoid it.

Should I change the function formula (I prefer not to), or should I revise the code, or change the solver?

Thank you in advance for any suggestion.

clear;
global G
tic
n=40;
A = ones(n) - eye(n); 
G = graph(A);

num_samples =1;
p.G = A
dim_G = size(p.G);
p.K = 1;
p.N = dim_G(1); 

nIters = 2000;
tBegin = 0;
tEnd = 100;

% initial condition
Init = -pi + 2*pi.*rand(p.N,num_samples);
dim_thetaInit = size(Init);
Init = Init(:,1);

[t,sol] = rk4(@fe,tBegin,tEnd,Init,nIters,p);

function td = fe(t,theta,p)
[theta_i,theta_j] = meshgrid(theta);
adj_matrix = p.G;
a = 4;
b = 1;
td= sum(adj_matrix'.*( ( b*cos(theta_i-theta_j) * (2*a^2*cos(theta_i-theta_j)*sin(theta_i-theta_j)-2*b^2*cos(theta_i-theta_j)*sin(theta_i-theta_j)) ) / ( 2*(b^2*(cos(theta_i-theta_j))^2+a^2*(sin(theta_i-theta_j))^2)^(3/2) ) + (b*sin(theta_i-theta_j))/(  sqrt( b^2*cos(theta_i-theta_j)^2 + a^2*sin(theta_i-theta_j)^2 )   )  ),1)';
end

function [T,Y] = rk4(f,a,b,ya,m,p)
%---------------------------------------------------------------------------
% RK4   Runge-Kutta solution for y' = f(t,y) with y(a) = ya .
% Sample call
%   [T,Y] = rk4('f',a,b,ya,m)
% Inputs
%   f    name of the function
%   a    left  endpoint of [a,b]
%   b    right endpoint of [a,b]
%   ya   initial value
%   m    number of steps
%   p    Kuramoto function arguments
% Return
%   T    solution: vector of abscissas
%   Y    solution: vector of ordinates
%
% NUMERICAL METHODS: MATLAB Programs, (c) John H . Mathews 1995
% To accompany the text:
% NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
% Prentice Hall, Englewood Cliffs, New Jersey, 07632, U . S . A .
% Prentice Hall, Inc .; USA, Canada, Mexico ISBN 0-13-624990-6
% Prentice Hall, International Editions:   ISBN 0-13-625047-5
% This free software is compliments of the author .
% E-mail address:      in % "mathews@fullerton.edu"
%
% Algorithm 9.4 (Runge-Kutta Method of Order 4) .
% Section   9.5, Runge-Kutta Methods, Page 460
%---------------------------------------------------------------------------
    h = (b - a)/m;
    T = zeros(1,m+1);
    size(ya,1)
    Y = zeros(size(ya,1),m+1);
    T(1) = a;
    Y(:,1) = ya;
    for j=1:m
      tj = T(j);
      yj = Y(:,j);
      k1 = h*feval(f,tj,yj,p);
      k2 = h*feval(f,tj+h/2,yj+k1/2,p);
      k3 = h*feval(f,tj+h/2,yj+k2/2,p);
      k4 = h*feval(f,tj+h,yj+k3,p);
      Y(:,j+1) = yj + (k1 + 2*k2 + 2*k3 + k4)/6;
      T(j+1) = a + h*j;
    end
end