MATLAB: How to write general function of 4th Order Runge-Kutta Method

errorindex out of boundsrunge kutta

I am trying to develop a Matlab function for the 4th Order Runge-Kutta Method. It needs to be able to work with any function for given initial conditions, step size, etc. and then plot the results afterwards. Here is the code that I have so far. There are no error in the function itself. However, when I try to use it, I get a couple of error messages that I have also shown below the mfile. I am not sure what the error messages mean or how I would go about correcting them. Any help would be much appreciated!
FUNCTION:
function [tp,yp] = rk4(dydt,tspan,y0,h)
if nargin < 4, error('at least 4 input arguments required'), end
if any(diff(tspan)<=0), error('tspan not ascending order'), end
n = length(tspan);
ti = tspan(1); tf = tspan(n);
tp = ti:h:tf;
yp = zeros(ti,length(tspan));
for n = ti:tf
    tp = ti;
    k1 = dydt(tp(n),y0(n));
    k2 = dydt(tp(n) + (1/2)*h, y0(n) + (1/2)*k1*h);
    k3 = dydt(tp(n) + (1/2)*h, y0(n) + (1/2)*k2*h);
    k4 = dydt(tp(n) + h, y0(n) + k3*h);
    yp = y0(n) + ((1/6)*(k1 + (2*k2) + (2*k3) + k4)*h);
    ti = ti + h;
    y0 = yp;
end
plot(tp,yp)
COMMAND WINDOW:
>> [tp,yp] = rk4(@(CA) ((1/5)*(20-CA))-(0.12*CA),15,0,1)
Attempted to access tp(15); index out of bounds because
numel(tp)=1.
Error in rk4 (line 13)
    k1 = dydt(tp(n),y0(n));

Best Answer

I didn’t run your code, but see if changing the initial ‘tp’ assignment to:
tv = ti:h:tf;
and then in your loop, changing:
tp = ti;
to:
tp = tv(n);
and changing all the subsequent references to ‘tp(n)’ to simply ‘tp’.
Also, you will need change this line (adding a subscript to ‘yp’) to save ‘yp’ at each iteration:
yp(n) = y0(n) + ((1/6)*(k1 + (2*k2) + (2*k3) + k4)*h);
See if that improves things. You may have to make some other changes as well to accommodate these, and there may be other problems in your code, but this should keep you progressing towards a successful conclusion.
Related Question