Why is this singular?

I'm trying to solve a system of differential equations. After a ridiculous amount of attempts I think i've finally managed to use ode15i correctly but it won't solve as the matrix is singular and I can't figure out why. I'd greatly appreacite it if someone could take a look at this and help me sort this out.

clc
n_total = 3.2203e21;
r = 2;  % [cm]
h = 30; % [cm]
Volume = pi*(r.^2)*h; % [cm^3]
e = 1.602e-19;% [C]
n_total = 3.22E21; % [1/cm^3], from P = nkT    

tspan = [0 1E-3];  
x0 = [1.5; 10E9; n_total; 0];
x01 = [0; 0 ;0 ;0];
%[X0,X01] = decic(@odefun,0,x0,[1 1 0 0],x01,[]);
[t,y] = ode15i(@odefun, tspan, x0, [0; 0; 0; 0])
nAr = y(:,3);
plot(t, nAr);

function res = odefun(t, x, x1) 
    r = 2;  % [cm]
    h = 30; % [cm]
    Volume = pi*(r.^2)*h; % [cm^3]
    e = 1.602e-19;% [C]
    kb = 1.38E-23;  % [J/K]
    Tg = 300;     % [K]
    Ti = Tg;
    PAr = 13.3322;   % [Pa], given as 0.1m[Torr]
    n_total = PAr/ (kb*Tg); % [1/cm^3], from P = nkT    
    A = r/2.405;
    uio = 1.53;    % [cm^2/ V-s]
    Patm = 101325; % [Pa], given as 1 atm
    natm = Patm/ (kb*Tg);
    Di = uio*(natm/n_total)*((kb*Ti)/e);   % [cm^2/ V-s],   n= n(total)   
    E1 = 11.6;   % [eV]
    E2 = 16;     % [eV] 

      Power_den = 100/Volume;

      res = zeros(4,1);
      k3 = 10000;
      k1 = ((2.5E-9).*(x(1).^0.74).*exp(-11.6./ x(1)));  %%[cm^3/s]
      k2 = ((2.3E-8).*(x(1).^0.68).*exp(-16./x(1)));     %%[cm^3/s]
      k4 = ((Di.*(1 + x(1)./ Ti)) ./A.^2);                                  %[1/s]

      Rate_sum = k1.*x(2).*x(3)*E1 + k2.*x(2).*x(3)*E2 + k3.*x(4) + k4.*x(2);

      res(1) = x1(1) - (2./ (3*x(2)*kb))*(Power_den - Rate_sum) + x1(2).*(x(1)./x(2));
      res(2) = x1(2) - x(2).*(k2.*x(3) - k4);
      res(3) = x1(3) + k1.*x(2).*x(3) + k2.*x(2).*x(3) - k3.*x(4) - k4.*x(2);
      res(4) = x1(4) - k1.*x(2).*x(3) + k3.*x(4);

end

Answer 1

[Torsten]

You return NaN values to ode15i within your vector "res". So you should check your inputs.

Best wishes

Torsten.