Saving vectors in a loop, wont save at time=0

Hi Guys,

I am currently trying to write code to solve ode's through the runge-kutta method. My code works fine for solving the odes, however i can't save the vector output from the ode's for the first loop. (iteration 0, time 0).

for t=0:h:tf

dx1=succinic_acid_production(t,x,..loads of other variables); K1=h*dx1; dx2=succinic_acid_production(t+h/2,x+K1/2,..loads of other variables); K2=h*dx2; dx3=succinic_acid_production(t+h/2,x+K2/2,k..loads of other variables); K3=h*dx3; dx4=succinic_acid_production(t+h,x+K3,..loads of other variables); K4=h*dx4;

    t=t+h;
    matrix(t,:)=x;
    x=x+(K1+2*K2+2*K3+K4)/6;

end

The solver works fine and gets me the correct curve when i plot(matrix), as shwon below, however the curves start at t=1. And i cant move the matrix(t,:)=x; section about t=t+h as otherwise i get the error message "Subscript indices must either be real positive integers or logicals."

I have been battling with this for far too long with no progress, any ideas?

Answer 1

[Mischa Kim]

Dan, I recommend to save the inital value of the DE before you enter the loop. Also because indexing in MATLAB starts at "1", which corresponds to t = 0. Something like:

 t = 0:h:tf;           % create time vector
 matrix(1,:) = x;      % save initial value

 for ii = 2:numel(t)   % use iteration steps as loop index rather than time

   dx1 = succinic_acid_production(t(ii),x,..loads of other variables); 
   K1  = h*dx1; 
   ...
   matrix(ii,:) = x;
   ...   
 end