title: How to solve an ODE 4th order with matlab ode23s?
body: <p>actually I am trying to solve the following ODE by using ode23s. Unfortunately I don't really know how to do it.</p>
  
  <p>The ode is given as: </p>
  
  <pre><code>bo(t)*y + b1(t)*y' + b2(t)*y'' + b3(t)*y''' + b4(t)*y''' = ao(t)*cst 
  </code></pre>
  
  <p>(Those derivatives should be marked with dots 'cause, they're all time-derivatives, eg. y' = dy/dt, y'' = d2y/dt2 etc. The parameters bo, b1, ..., b4 and ao are given as matlab-functions. 'cst' is a constant.</p>
  
  <p>First, I transformed this equation into a system of ODEs, what leads me to the following function:</p>
  
  <pre><code>function [dy] = odetest(y,cst,ao,bo,b1,b2,b3,b4)
    dy = [ y(2); ...
           y(3); ...
           y(4); ...
           (-b1*y(2) + ao*cst - b2*y(3) - b3*y(4) + bo*y(1))*(1/(b4)) ];
  end
  </code></pre>
  
  <p>So ... to simplify the procedure (and as a first approximation) I want to hold the time-dependent parameters ao, bo, b1, ... b4 constant for each time step.</p>
  
  <p>And now the problem starts: I don't know how to put this function in ode23s that it could be solved for an arbitrary timestep with given initial-conditions. Maybe I made a kind of mistake in the transformation of the time-derivatives ...</p>
  
  <p>I guessed that the solution has to be something like, which unfortunately doesn't work</p>
  
  <pre><code>% file where the Ode has to be solved
  % Loop over all time-steps
  
  for i=1:number_timesteps
  
     time0 = i-1;
     time  = i;
  
     ao = function_ao(i); % function of parameter a, which is now constant in each timestep
  
     bo = function_bo(i);
     b1 = function_b1(i);
     b2 = function_b2(i);
     b3 = function_b3(i);
     b4 = function_b4(i);
  
     y = value1;   % initial conditions for y
     dy = value2;
     ddy = value3;
     dddy = value4;
  
     cst = value5; % value of cst in the actual timestep
  
     [T, Y] = ode23s(@odetest, [time0 time1] [y; dy; ddy; dddy], ...
     cst, ao, bo, b1, b2, b3, b4 ); % the arguments are the inital conditions, parameters ao, bo, ...
  
  end
  </code></pre>
  
  <p>But ... I don't really understand what I did wrong. And I don't know either how to solve the problem if the parameters ao, bo, b1, ..., b4 change in one time-step, so that ode23s has to use their functions somehow.</p>
  
  <p>It would be great if someone can help me. Thank you in advance!</p>
  
  <p><strong>Edit - Implementation</strong></p>
  
  <p>I now modified my code as follows, but unfortunately get an error "Error using odetest (line 5) Not enough input arguments." I don't know why :/</p>
  
  <p>odetest.m:</p>
  
  <pre><code>    function [ dy ] = odetest( ttime, y, cst, ff1, ff2, ff3, ff4, ...
                    df1, df2, df3, df4, dd2, dd31, dd32, dd33, ...
                    dd34, ttk, tt0, EE1)
  
      EE0 = fe0(ff1, ff2, ff3, ff4, ttk, tt0, ttime);
      e1 = fe1(df1, dd2, dd31, ttk, tt0, ttime);
      e2 = fe2(df2, dd2, dd32, ttk, tt0, ttime);
      e3 = fe3(df3, dd2, dd33, ttk, tt0, ttime);
      e4 = fe4(df4, dd2, dd34, ttk, tt0, ttime);
  
      a0 = fa0(EE0, e1, e2, e3, e4);    
      b0 = fb0(e1, e2, e3, e4);
      b1 = fb1(EE0, EE1, EE1, EE1, EE1, e1, e2, e3, e4);
      b2 = fb2(EE0, EE1, EE1, EE1, EE1, e1, e2, e3, e4);
      b3 = fb3(EE0, EE1, EE1, EE1, EE1, e1, e2, e3, e4);
      b4 = fb4(EE0, EE1, EE1, EE1, EE1);
  
      dy = [ y(2); ...
             y(3); ...
             y(4); ...
             (-bb_1*y(2) + aa_0*cst - bb_2*y(3) - bb_3*y(4) + bb_0*y(1))*(1/(bb_4))];
  
  end
  </code></pre>
  
  <p>The file itself, where I call this function is</p>
  
  <pre><code>% something before (set constant parameters val_i and define initial
  % conditions yy_t0 for first step ...
  for i=1:n
  
     dbl_t0 = t(i-1,1);
     dbl_t  = t(i,1);
  
     cst = act(i,1) - act(i-1);
  
  
     [T, Y] = ode23s(@odetest, [dbl_t0 dbl_t], [yy_t0 ...
              dyy_t0 ddy_t0 dddy_t0],cst, ...
              val_f1, val_f2, val_f3, val_f4, val_d1, val_d2, val_d3, val_d4, val_dd2, ...
              val_dk1, val_dk2, val_dk3, val_dk4, val_tk, val_t0, val_E1);
  
     yy_t = Y(1);
     dyy_t = Y(2);
     ddyy_t = Y(3);
     dddyy_t = Y(4);
  
     yy_t0 = yy_t;
     dyy_t0 = dyy_t;
     ddyy_t0 = ddyy_t;
     dddyy_t0 = dddyy_t;
  
  end
  </code></pre>
  
  <p>... does anyone sees a mistake in here?</p>
  
  <p><strong>Last Edit</strong></p>
  
  <p>I found the mistake. - I have to call ode23s in the following way:</p>
  
  <pre><code>[T, Y] = ode23s(@(ttime,y)odetest(ttime, y, cst, ...
              f_1, f_2, f_3, f_4, d_11, d_12, d_13, d_14, d_2, d_31, ...
              d_32, d_33, d_34, t_k, t_0, E_1), [dbl_t0 dbl_t], ...
              [yy_t0 dyy_t0 ddyy_t0 dddyy_t0]);
  </code></pre>
  
  <p>... now it works.</p>