function J = SIR_fun(q,Y0,t_data,I_data)
  
  ode_options = odeset('RelTol',1e-6);
  
  [t,Y] = ode45(@SIR_rhs,t_data,Y0,ode_options,q);
  I = Y(:,2);
  Res = I - I_data;
  J = Res'*Res;