Solving 1st ODE, while cycling a solution into 2nd ODE and using that solution in 1st ODE

I'm having an issue making this work with a set of reaction equations, which I've set the code up for:

function dcdt = Rxn( t, C )
  global Tsmat
  % 
  % Runs at a set time point
  % f [=] mol/s , w [=] mol-site
  k = [ 4.21E13, 1.02E15, 4.68E10, 4.21E13, 1.19E9, 6.2E10, 6.1E9, 1.8E5, 1.23E5, 9.87E-1, 4.61E1, 5.05, 2.27E3, 2.25E3, 5.05 ]; %reaction coeffcients
  E = [ 111450, 129530, 165160, 111450, 52374, 90063, 69237, 56720, 81920, 5296, 25101, 21768, 39070, 39680, 21768 ]; % activation energies
  kG = [ 7.33, 2.57E2, 1.8E-4, 3.14E6 ]; % parameters f or solving inhibition factor G
  EG = [ -485, 166, -10163, 3685 ];
  Rg = 8.314; % [=] j/mol/K
  z = 1; % [=] K
  P = 1000000; % [=] Pascals J/m^3
  Ts = Tsmat ( z, t );
  Csi=P/Rg/Ts;
  Gwgs = -4.1034E4 + 44.19 * Ts - 5.553E-3 * Ts ^ 2;
  Kwgs = exp ( -Gwgs / ( Rg * Ts ) );
  gval = [ kG; EG ];
  K = zeros ( 4 );
  theta=1;
  for i = 1:4
    K (i) = gval ( 1, i ) * exp ( -gval( 2, i ) / Ts );
  end
 G = ( ( 1 + K ( 1 ) * C ( 1 ) + K ( 2 ) * C ( 3 ) ).^2 ) * ( 1 + K ( 3 ) * C ( 1 ).^2 * C ( 3 ).^2 ) * ( 1 + K ( 4 ) * C( 8 ) );
 %From C(1) to C(11) [Co, O2, C3H6, C3H8, H2, Co2, H2O, NO, N2, Ce2O3, Ceo2]
 R1 = k ( 1 ) * C ( 1 ) * C ( 2 ) * exp ( -E( 1 ) / ( Rg * Ts ) ) / G; %[=] mol/mol-site/s
 R2 = k ( 2 ) * C ( 3 ) * C ( 2 ) * exp ( -E( 2 ) / ( Rg * Ts ) ) / G;
 R3 = k ( 3 ) * C ( 4 ) * C ( 2 ) * exp ( -E( 3 ) / ( Rg * Ts ) ) / G;
 R4 = k ( 4 ) * C ( 5 ) * C ( 2 ) * exp ( -E ( 4 ) / ( Rg * Ts ) ) / G;
 R5 = k ( 5 ) * C ( 1 ) * C ( 8 ) * exp ( -E ( 5 ) / ( Rg * Ts ) ) / G;
 R6 = k ( 6 ) * C ( 3 ) * C ( 8 ) * exp ( -E ( 6 ) / ( Rg * Ts ) ) / G;
 R7 = k ( 7 ) * C ( 5 ) * C ( 8 ) * exp ( -E ( 7 ) / ( Rg * Ts ) ) / G;
 R8 = k ( 8 ) * ( C ( 1 ) * C ( 7 ) - C ( 5 ) * C ( 6 ) / Kwgs ) * exp ( -E ( 8 ) / ( Rg * Ts ) ) / G;
 R9 = k ( 9 ) * C ( 3 ) * C ( 7 ) * exp ( E ( 9 ) / ( Rg * Ts ) ) / G;
 R10 = k ( 10 ) * C ( 2 ) * exp ( -E ( 10 ) / ( Rg * Ts ) ) * ( 1 - theta );
 R11 = k ( 11 ) * C ( 8 ) * exp ( -E ( 11 ) / ( Rg * Ts ) ) * ( 1 - theta );
 R12 = k ( 12 ) * C ( 1 ) * exp ( -E ( 12 ) / ( Rg * Ts ) ) * theta;
 R13 = k ( 13 ) * C ( 3 ) * exp ( -E ( 13 ) / ( Rg * Ts ) ) * theta;
 R14 = k ( 14 ) * C ( 4 ) * exp ( -E ( 14 ) / ( Rg * Ts ) ) * theta;
 R15 = k ( 15 ) * C ( 5 ) * exp ( -E ( 15 ) / ( Rg * Ts ) ) * theta;
% Given in R [=] mol/mol-site/sec, multipy by aj to get mol/m^3
 dcdt ( 1, 1 ) = -R1 - R5 - R8 + 3 * R9 - R12 + 3 * R13 + 3 * R14; %CO [=] mol/mol-site/s
 dcdt ( 2, 1 ) = -0.5 * R1 - 4.5 * R2 - 5 * R3 - 0.5 * R4 - R10; %O2
 dcdt ( 3, 1 ) = -R2 - R6 - R9 - R13; %C3H6
 dcdt ( 4, 1 ) = -R3 - R14; %C3H8
 dcdt ( 5, 1 ) = -R4 - R7 + R8 + 6 * R9 - R15; %H2
 dcdt ( 6, 1 ) = R1 + 3 * R2 + 3 * R3 + R5 + 3 * R6 + R8 + R12; %C(2)
 dcdt ( 7, 1 ) = 3 * R2 + 4 * R3 + R4 + 3 * R6 + R7 - R8 - 3 * R9 + 3 * R13 + 4 * R14 + R15; %H2O
 dcdt ( 8, 1 ) = -R5 - 9 * R6 - R7 - R11; %NO
 dcdt ( 9, 1 ) = 0.5 * R5 + 4.5 * R6 + 0.5 * R7 + 0.5 * R11; %N2
 dcdt ( 10, 1 ) = -2 * R10 - R11 + R12 + 6 * R13 + 7 * R14 + R15; %Ce2O3
 dcdt ( 11, 1 ) = 4 * R10 + 2 * R11 - 2 * R12 - 12 * R13 - 14 * R14 - 2 * R15;%CeO2
end

I'm trying to work in this equation

Any ideas? I've tried simplifying the above equation using finite difference and then solving the above code simultaneously with this code,

function theta_new = dthetadt (theta,tspan,c1,c2,c3,c4,c5,c6, i)

    t_num = 201;
    t_min = 0.0;
    t_max = 80.0;
    dt = ( t_max - t_min ) / ( t_num - 1 );
    if i == 2
    theta_new = zeros(t_num,1);
    end

    theta_new (i) = ((4*c1+2*c2)-(2*c3+12*c4+14*c5+2*c6))*dt+theta(i-1);

    return
  end

with the edited first code:

function dcdt = Rxn( t, C )
  global Tsmat
  % 
  % Runs at a set time point
  % f [=] mol/s , w [=] mol-site
  k = [ 4.21E13, 1.02E15, 4.68E10, 4.21E13, 1.19E9, 6.2E10, 6.1E9, 1.8E5, 1.23E5, 9.87E-1, 4.61E1, 5.05, 2.27E3, 2.25E3, 5.05 ]; %reaction coeffcients
  E = [ 111450, 129530, 165160, 111450, 52374, 90063, 69237, 56720, 81920, 5296, 25101, 21768, 39070, 39680, 21768 ]; % activation energies
  kG = [ 7.33, 2.57E2, 1.8E-4, 3.14E6 ]; % parameters f or solving inhibition factor G
  EG = [ -485, 166, -10163, 3685 ];
  Rg = 8.314; % [=] j/mol/K
  z = 1; % [=] K
  P = 1000000; % [=] Pascals J/m^3
  Ts = Tsmat ( z, t );
  Csi=P/Rg/Ts;
  Gwgs = -4.1034E4 + 44.19 * Ts - 5.553E-3 * Ts ^ 2;
  Kwgs = exp ( -Gwgs / ( Rg * Ts ) );
  gval = [ kG; EG ];
  K = zeros ( 4 );
  for t == 1
      theta = zeros(201,1);
  end
  theta(1)=1;
  for i = 1:4
    K (i) = gval ( 1, i ) * exp ( -gval( 2, i ) / Ts );
  end
 G = ( ( 1 + K ( 1 ) * C ( 1 ) + K ( 2 ) * C ( 3 ) ).^2 ) * ( 1 + K ( 3 ) * C ( 1 ).^2 * C ( 3 ).^2 ) * ( 1 + K ( 4 ) * C( 8 ) );
 %From C(1) to C(11) [Co, O2, C3H6, C3H8, H2, Co2, H2O, NO, N2, Ce2O3, Ceo2]
 R1 = k ( 1 ) * C ( 1 ) * C ( 2 ) * exp ( -E( 1 ) / ( Rg * Ts ) ) / G; %[=] mol/mol-site/s
 R2 = k ( 2 ) * C ( 3 ) * C ( 2 ) * exp ( -E( 2 ) / ( Rg * Ts ) ) / G;
 R3 = k ( 3 ) * C ( 4 ) * C ( 2 ) * exp ( -E( 3 ) / ( Rg * Ts ) ) / G;
 R4 = k ( 4 ) * C ( 5 ) * C ( 2 ) * exp ( -E ( 4 ) / ( Rg * Ts ) ) / G;
 R5 = k ( 5 ) * C ( 1 ) * C ( 8 ) * exp ( -E ( 5 ) / ( Rg * Ts ) ) / G;
 R6 = k ( 6 ) * C ( 3 ) * C ( 8 ) * exp ( -E ( 6 ) / ( Rg * Ts ) ) / G;
 R7 = k ( 7 ) * C ( 5 ) * C ( 8 ) * exp ( -E ( 7 ) / ( Rg * Ts ) ) / G;
 R8 = k ( 8 ) * ( C ( 1 ) * C ( 7 ) - C ( 5 ) * C ( 6 ) / Kwgs ) * exp ( -E ( 8 ) / ( Rg * Ts ) ) / G;
 R9 = k ( 9 ) * C ( 3 ) * C ( 7 ) * exp ( E ( 9 ) / ( Rg * Ts ) ) / G;
 R10 = k ( 10 ) * C ( 2 ) * exp ( -E ( 10 ) / ( Rg * Ts ) ) * ( 1 - theta(t) );
 R11 = k ( 11 ) * C ( 8 ) * exp ( -E ( 11 ) / ( Rg * Ts ) ) * ( 1 - theta(t) );
 R12 = k ( 12 ) * C ( 1 ) * exp ( -E ( 12 ) / ( Rg * Ts ) ) * theta(t);
 R13 = k ( 13 ) * C ( 3 ) * exp ( -E ( 13 ) / ( Rg * Ts ) ) * theta(t);
 R14 = k ( 14 ) * C ( 4 ) * exp ( -E ( 14 ) / ( Rg * Ts ) ) * theta(t);
 R15 = k ( 15 ) * C ( 5 ) * exp ( -E ( 15 ) / ( Rg * Ts ) ) * theta(t);
% Given in R [=] mol/mol-site/sec, multipy by aj to get mol/m^3
 for i = 2:201
    theta(i) = dthetadt(theta, tspan, c1, c2, c3, c4, c5, c6, i);
 end
 dcdt ( 1, 1 ) = -R1 - R5 - R8 + 3 * R9 - R12 + 3 * R13 + 3 * R14; %CO [=] mol/mol-site/s
 dcdt ( 2, 1 ) = -0.5 * R1 - 4.5 * R2 - 5 * R3 - 0.5 * R4 - R10; %O2
 dcdt ( 3, 1 ) = -R2 - R6 - R9 - R13; %C3H6
 dcdt ( 4, 1 ) = -R3 - R14; %C3H8
 dcdt ( 5, 1 ) = -R4 - R7 + R8 + 6 * R9 - R15; %H2
 dcdt ( 6, 1 ) = R1 + 3 * R2 + 3 * R3 + R5 + 3 * R6 + R8 + R12; %C(2)
 dcdt ( 7, 1 ) = 3 * R2 + 4 * R3 + R4 + 3 * R6 + R7 - R8 - 3 * R9 + 3 * R13 + 4 * R14 + R15; %H2O
 dcdt ( 8, 1 ) = -R5 - 9 * R6 - R7 - R11; %NO
 dcdt ( 9, 1 ) = 0.5 * R5 + 4.5 * R6 + 0.5 * R7 + 0.5 * R11; %N2
 dcdt ( 10, 1 ) = -2 * R10 - R11 + R12 + 6 * R13 + 7 * R14 + R15; %Ce2O3
 dcdt ( 11, 1 ) = 4 * R10 + 2 * R11 - 2 * R12 - 12 * R13 - 14 * R14 - 2 * R15;%CeO2
end

however I got the error:

Error using ode1 (line 40)
Unable to evaluate the ODEFUN at t0,y0. Error: File: Rxn.m Line: 19 Column: 9
Unexpected MATLAB operator.

Error in runrxn (line 15)
C=ode1(@Rxn,t,[10,10,3,0,0,0,10,5,8,0,0]);

Any advice would be greatly appreciated!

[Torsten]

Replace "for t==1" by "if t==1".

But you will run into difficulties because the solver usually won't hit t=1 exactly.

Best wishes

Torsten.

Answer 1

[Jan]

See Torsten's comment: I'm not sure if this is the only problem, but this should cause problems:

for t == 1   % See error message: == is the elementwise comparison
             % A loop would be: for t = 1
             % But a loop over 1 element is not meaningful here.
  theta = zeros(201,1);
end
theta(1) = 1;

Now theta is [1] except if the integrator hits t=1.0 by accident, that theta is [1, 0,0,0...]. This cannot be wanted.

theta(t) cannot work also, because t is not a positive integer in general.

[Cpan]

In that case, when using ode45 is there any way to introduce a variable within the function that I'm trying to solve and replace it with each ode45 iteration? Since within my code I need to use values within the function that are solved within ode45, but I do have an initial condition for theta.

[Jan]

This sounds, like it will cause discontinuities in the function to be integrated. Then ODE45 cannot handle the function correctly.

What does "in each ode45 iteration" mean? A step can be rejected by the step size control, but is this "an iteration" or not?

I'm not sure, what you exactly want. Can you express it mathematically?

[Cpan]

I realized that I could just pass the initial condition as C(12) and then ode45 would take care of the rest. Huge brain fart on my part. Thank you for taking a look at my question.

[Walter Roberson]

Cpan, see

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html