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

I am solving a system of ODE to model electrospinning, which could be described as n Viscoelastically connected point charge

    y0=[r(LeadingBeadidx:LastBeadidx+1,:),v(LeadingBeadidx:LastBeadidx+1,:),...
        sigmau(LeadingBeadidx:LastBeadidx+1)];
    y0=reshape(y0,[],1);
    %ODE45 expects a column vector
    [t,YSol]=ode45(@(t,y) MasterODE(t,y,m(LeadingBeadidx:LastBeadidx+1),...
        e(LeadingBeadidx:LastBeadidx+1),alpha,E,G,mu,InitSegV),...
        [0,dt],y0,ODEoptions);
    YSol=reshape(YSol(t==dt,:),[],7);
function dVardt= MasterODE(t,Var,m,e,alpha,E,G,mu,InitSegV) %#ok<INUSL>
%Mater ODE that combines motion(newton law) and sigma(maxwell materials for
%all beads

GPUCompute=numel(Var)>=700 && gpuDeviceCount>=1;
if GPUCompute
    try
        %transfer larger data to GPU, if present
        Var=gpuArray(Var);
        m=gpuArray(m);
        e=gpuArray(e);
    catch
    end
end
Var=reshape(Var,[],7);
r=Var(:,1:3);
v=Var(:,4:6);
sigmau=Var(:,7);

%========================================================================
%calculate required parameters from r and v for force calculation
%========================================================================

f=Viscoelastic(au,sigmau,dispu,lu)-Viscoelastic(ad,sigmad,dispd,ld);
%Calculate f, starting with viscoelastic forces

f=f+SurfaceTension(au,ad,r,dispcross,dispu,dispd,alpha);
%adding surface tension

f=f+ExternalField(e,E);
%Adding external electric field

f=f+InterElectric(e,r);
%Adding inter beads electrostatic forces
%========================================================================

drdt=v;
%Derivative of position is velocity

dvdt=f./m;
%f=ma,a=dvdt

dsigmaudt=(G./lu).*dludt-(G./mu).*sigmau;
%Derivatives of sigma

dVardt=[drdt,dvdt,dsigmaudt];
dVardt(isnan(dVardt))=0;
%Assemble derivatives

dVardt=reshape(dVardt,[],1);%ODE45 expects column vector

if GPUCompute
    dVardt=gather(dVardt);% Return array to workspace
end
end

where r and v are n-by-3 vectors for the position and velocity of n points, sigma is a n-by-1 vector of stress between points.

This system could be solved by ode45 reasonably fast when the n is small. However when n is large (like >1000) it takes forever, even for a small time span. Feeling that the system may become stiff, I was tempted to try some stiff solvers like ode15s, but got greeted by the following error:

...
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. 
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. 
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. 
Warning: Failure at t=0.000000e+00.  Unable to meet integration tolerances without reducing the step size below the smallest
value allowed (7.905050e-323) at time t. 

So my questions are:

1.What is the cause of this error? How do I solve it?

2.I could supply part of the jacobian, or at least a jPattern. Will this help?

3. Does this error means that this is not a stiff problem and not a good fit for stiff solvers? How does one determine when the system is stiff?

Thanks!

Answer 1

[Yi-xiao Liu]

OK, I am going to answer my own question here.

After providing Jpattern, the ode15s works normally, but does not display any performance advantages over ode45. But at least Jpattern makes it work.