Why am i getting blank plots

I cant seem to get a plot and i cant figure out why. This is my code function

clear clc close all
% Exothermic Reactor solution file
% Define the System Inputs 
SI.kref = 0.5; % Arrhenius constant (/s) 
SI.E = 20e3; % activation energy (J/mol) 
SI.R = 8.314; % ideal gas constant (J/mol/K) 
SI.M1 = 200; % mass of material in CSTR 1 (kg) 
SI.M2 = 300; % mass of material in CSTR 2 (kg) 
SI.cp = 4180; % Specific heat capcity in the reactor (J/kg/oC) 
SI.theta_Hi = 100; % temperature of jacket fluid in (oC) 
SI.theta_F = 50; % initial temperature of Feed (oC) 
SI.theta_ref = 20; % reference temperature (oC) 
SI.F = 5; % flowrate between CSTRs (kg/s) 
SI.q = 0.5; % flow rate of heating fluid in jackets (kg/s) 
SI.cph = 4000; % Specific heat capcity of heating fluid (J/kg/oC) 
SI.CF = 1; % Feed Concentration to CSTR 1 (mol/kg) 
SI.U1 = 300; % overall heat transfer coefficient to CSTR 1 (W/m^2/oC) 
SI.U2 = 300; % overall heat transfer coefficient to CSTR 2 (W/m^2/oC) 
SI.A1 = 3; % jacket 1 internal heat transfer area (m^2) 
SI.A2 = 3; % jacket 2 internal heat transfer area (m^2)
% Define Initial Values for ODES 
theta_1i = 20; % (oC) 
theta_2i = 20; % (oC) 
C1_i = 10; % (mol/kg) 
C2_i = 10; % (mol/kg) 
D_i = [theta_1i; theta_2i; C1_i; C2_i]; % vector of initial conditions
% Define the time for the simulation 
t_end = 3*3600; % s 
t_sim = 0:60:t_end; % s
% Solve the ODEs using ode45 
[t, D] = ode45(@(t,D) ExothermicODEs(t,D,SI),t_sim, D_i);
theta_1 = D(:,1) 

theta_1 = 181×1

20 NaN NaN NaN NaN NaN NaN NaN NaN NaN

theta_2 = D(:,2) 

theta_2 = 181×1

20 NaN NaN NaN NaN NaN NaN NaN NaN NaN

C1 = D(:,3) 

C1 = 181×1

10 NaN NaN NaN NaN NaN NaN NaN NaN NaN

C2 = D(:,4)

C2 = 181×1

10 NaN NaN NaN NaN NaN NaN NaN NaN NaN

figure % Graph of temperature vs time for CSTR 1 
subplot(2,2,1), plot(t, theta_1, 'b-'); 
xlabel('Time (s)','FontSize', 16); 
ylabel('Temperature (\circC)', 'FontSize', 16);
% Graph of temperature vs time for CSTR 2 
subplot(2,2,2), plot(t, theta_2, 'b-'); 
xlabel('Time (s)','FontSize', 16); 
ylabel('Temperature (\circC)', 'FontSize', 16);
% Graph of Concentration vs time for CSTR 1 
subplot(2,2,3), plot(t, C1, 'g-'); 
xlabel('Time (s)','FontSize', 16); 
ylabel('Concentration (mol/kg)', 'FontSize', 16);
% Graph of Concentration vs time for CSTR 2 
subplot(2,2,4), plot(t, C2, 'g-'); 
xlabel('Time (s)','FontSize', 16); 
ylabel('Concentration (mol/kg)', 'FontSize', 16);

function dD_dt = ExothermicODEs(~,D,SI) % Problem 1 ODEs 
  theta_1 = D(1); % CSTR 1 temperature (oC) 
  theta_2 = D(2); % CSTR 2 temperature (oC) 
  C1 = D(3); % CSTR 1 concentration (mol/kg) 
  C2 = D(4); % CSTR 2 concentration (mol/kg)
  M1 = SI.M1; % mass of material in CSTR 1 (kg) 
  M2 = SI.M2; % mass of material in CSTR 2 (kg) 
  cp = SI.cp; % Specific heat capcity in the reactor (J/kg/oC) 
  theta_Hi = SI.theta_Hi; % temperature of jacket fluid in (oC) 
  theta_F = SI.theta_F; % initial temperature of Feed (oC) 
  E = SI.E; % activation energy (J/mol) 
  R = SI.R; % ideal gas constant (J/mol/K) 
  kref = SI.kref; % Arrhenius Constant (/S) 
  F = SI.F; % flowrate between CSTRs (kg/s) 
  q = SI.q; % flow rate of heating fluid in jackets (kg/s) 
  cph = SI.cph; % Specific heat capcity of heating fluid (J/kg/oC) 
  CF = SI.CF; % Feed Concentration to CSTR 1 (mol/kg) 
  U1 = SI.U1; % overall heat transfer coefficient to CSTR 1 (W/m^2/oC) 
  U2 = SI.U2; % overall heat transfer coefficient to CSTR 2 (W/m^2/oC) 
  A1 = SI.A1; % jacket 1 internal heat transfer area (m^2) 
  A2 = SI.A2; % jacket 2 internal heat transfer area (m^2) 
  theta_ref = SI.theta_ref; % referance temperature (oC)
  % algebraic eqautions 
  k1 = kref*exp((E/R)*((1/(theta_ref+273.15))-(1/(theta_1+273.15)))); % rate constant in CSTR 1 (/s) 
  k2 = kref*exp((E/R)*((1/(theta_ref+273.15))-(1/(theta_2+273.15)))); % rate constant in CSTR 2 (/s) 
  theta_Ho = theta_1 + (theta_Hi-theta_1)*exp((-U1*A1)/(q*cph)); % temperature of jacket fluid out (oC) 
  theta_Hm = theta_2 + (theta_Hi-theta_2)*exp((-U2*A2)/(q*cph)); % temperature of jacket between CSTR 1 and 2 (oC) 
  theta_lm1 = (theta_Hm-theta_Ho)/(log((theta_Hm-theta_1)/(theta_Ho-theta_1))); % Log mean temperature difference for CSTR 1 (oC) 
  theta_lm2 = (theta_Hm-theta_Ho)/(log((theta_Hm-theta_2)/(theta_Ho-theta_2))); % Log mean temperature difference for CSTR 2 (oC)
  % ODEs 
  dtheta_1_dt = ((F*(theta_F-theta_1))/M1) + ((U1*A1*theta_lm1)/(M1*cp)); % CSTR 1 temperature ODE 
  dtheta_2_dt = ((F*(theta_F-theta_2))/M2) + ((U2*A2*theta_lm2)/(M2*cp)); % CSTR 2 temperature ODE 
  dC1_dt = ((F*(CF-C2))/M1)-k1*C1; % CSTR 1 concentration ODE 
  dC2_dt = ((F*(C1-C2))/M2)-k2*C2; % CSTR 2 concentration ODE
  dD_dt = [dtheta_1_dt;dtheta_2_dt; dC1_dt; dC2_dt];
end

Answer 1

[Torsten]

The solution for all four solution variables is NaN (see above) - most probably caused by the use of "log" in "theta_lm1" and "theta_lm2".