title: Splitting the state vector when solving a series of odes using python odeint scipy
body: <p>I'm converting my code from Matlab to Python and stuck on how to split the state vector such that the result returns the two solution. I have a vector and a single value for the two initial conditions and I expect as the final result a matrix and a vector.</p>
  
  <p>I tried joining the initial conditions (y0 = [c_pt_0, x_0]) in the same manner as the solution (soln = [dfdt,dcdt]) (which is shown below in the code). I also tried a similar approach that is used in matlab, which is concatenating the the initial conditions to one single array and unpacking the results but I think the problem is in the dimensions. </p>
  
  <pre><code>#Basic imports
  import numpy as np
  import pylab
  import matplotlib. pyplot as plt
  import scipy
  from scipy.integrate import odeint
  import matplotlib.pyplot as plt
  
  # define parameters
  pi = 3.14159265 
  V_m = 9.09
  m_V__M_Pt = 1e6/195.084    
  rho = 21.45
  R0 = 10**(-8.19)
  k_d = 10**(-13)
  k_r = 10**(-5)
  S = 0.314   #distribution parameter
  M = 0.944   #distribution parameter
  
  ## geometry
  # Finite Volume Method with equidistant elements
  r_max = 30.1e-9                     #maximum value
  n = 301                             #number of elements of FVM
  dr = r_max/n                        #length of elements, equidistant
  ini_r = np.linspace(5e-10,r_max,n+1)    #boundaries of elements
  mid_r = ini_r[0:n]+dr/2             #center of elements
  
  ## initial conditions
  #initial distribution
  x0 = 1/(S*np.sqrt(2*pi)*mid_r*1e9)*np.exp((-(np.log(mid_r*1e9)-M)**2)/(2*S**2))
  c_pt_0 = 0
  y0 = [x0, c_pt_0]
  
  MN_0 = scipy.trapz(np.power(mid_r, 3)*x0,
                          x=mid_r)    # initial mass
  M_0 = 4/3*pi*rho*MN_0
  
  def f(y, t):
      r = y[0]
      c_pt = y[1]
  
      #materials balance
      drdt = V_m * k_r * c_pt * np.exp(-R0/ mid_r) - V_m * k_d * np.exp(R0/ mid_r)
      dmdt = 4*pi*rho*mid_r**2*drdt
      dMdt = np.trapz(r*dmdt, x=mid_r)
  
      dcdt = m_V__M_Pt*(-dMdt)/M_0
      dfdt = -(np.gradient(r*drdt, dr))
  
      soln = [dfdt, dcdt]
      return soln
  
  #------------------------------------------------------
  #define timespace
  time = np.linspace(0, 30, 500)
  
  #solve ode system
  sln_1 = odeint (f , y0 , time,
      rtol = 1e-3, atol = 1e-5)
  
  pylab.plot(mid_r, sln_1[1,:], color = 'r', marker = 'o')
  pylab.plot(mid_r, sln_1[-1,:], color = 'b', marker = 'o')
  plt.show()
  </code></pre>
  
  <p>Traceback:
  ValueError: setting an array element with a sequence.</p>
  
  <p>Any help is much appreciated.</p>
  
  <p><strong>EDIT: ADDED MATLAB CODE</strong></p>
  
  <p>Here is the MATLAB code that works that I want to convert to python where the state vector is split. I have three files (one main, the f function, and the parameters). Please excuse any <em>face palm</em> coding errors but I do appreciate any suggestions even for this.</p>
  
  <p>modified_model.m:</p>
  
  <pre><code>function modified_model
  
  % import parameters
  p = cycling_parameters;
  
  % initial conditions
  c_pt_0 = 0;
  y0 = [p.x0; c_pt_0];
  
  % call integrator 
  options_ODE=odeset('Stats','on', 'RelTol',1e-3,'AbsTol',1e-5);
  [~, y] = ode15s(@(t,y) f(t, y, p), p.time, y0, options_ODE);
  
  %% Post processing
  % split state vector
  r       = y(:,1:p.n);
  c_Pt = y(:,p.n+1);
  
  %% Plot results 
  
  figure
  hold on;
  plot(p.r_m, r(1,:));
  plot(p.r_m, r(end,:));
  xlabel({'size'},'FontSize',15)
  ylabel({'counts'},'FontSize',15)
  </code></pre>
  
  <p>f.m</p>
  
  <pre><code>function soln = f(~, y, p)
  
  %split state vector
  r = y(1:p.n);
  c_pt = y(p.n+1);
  
  % materials balance
  drdt = p.Vm_Pt.*p.k_rdp.*c_pt.*exp(-p.R0./p.r_m) - p.Vm_Pt.*p.k_dis.*exp(p.R0./p.r_m);
  dmdt = 4*pi*p.rho*p.r_m.^2.*drdt; 
  dMdt = trapz(p.r_m, r.*dmdt);
  dcdt = p.I_V*p.m_V__M_Pt*(-dMdt)/p.M_0;
  dfdt = - gradient(r.*drdt,p.dr);
  
  soln = [dfdt; dcdt];
  </code></pre>
  
  <p>and the parameters file: cycling_parameters.m</p>
  
  <pre><code>function p=cycling_parameters
  
  p.M = 195.084; 
  p.rho = 21.45;       
  p.time = linspace(0, 30, 500);
  p.m_V__M_Pt = 1e6/p.M;     
  
  p.Vm_Pt = 9.09;  
  
  p.R0_log = -8.1963; 
  p.k_dis_log = -13; 
  p.k_rdp_log = -11;
  
  p.R0 = 10^(p.R0_log);
  p.k_dis = 10^(p.k_dis_log);
  p.k_rdp = 10^(p.k_rdp_log);
  
  %%% geometry %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%% Finite Volume Method with equidistant elements 
  p.r_max = 10.1e-9;                    % [m] maximum radius of PRD
  
  p.n = 301;                          % number of elements of FVM
  
  p.dr = p.r_max/p.n;                 % [m] length of elements, equidistant
  p.r = linspace(5e-10,p.r_max,p.n+1)';   % [m] boundaries of elements
  p.r_m = p.r(1:p.n)+p.dr/2;          % [m] center of elements
  
  %log normal initial distribution
  
  S = 0.314; 
  M = 0.944; 
  p.x0 = 1./(S.*sqrt(2.*pi).*p.r_m*1e9).*exp((-(log(p.r_m*1e9)-M).^2)./(2.*S.^2));
  
  
  p.r_squared = p.r_m.^2; % [m^2] squares of the radius (center of elements)
  p.r_cubed   = p.r_m.^3; % [m^3] cubes of the radius (center of elements)
  p.MN_0 = trapz(p.r_m, p.r_cubed.*p.x0);   % Eq. 2.11 denominator
  p.M_0 = 4/3*pi*p.rho*p.MN_0;  
  p.I_V = 1; %ionomer volume fraction in the catalyst layer
  </code></pre>