Implementation of 4-step Runge-Kutta

I am trying to use a runge-kutta function that I wrote previously on a new example.

function [y, t] = rungekutta(f, a, b, y0, h)

t = a:h:b;

y = zeros(size(t));

y(1) = y0;

for i = 1:length(t)-1

k1 = f(t(i), y(i));

k2 = f(t(i) + h/2, y(i) + (h/2)*k1);

k3 = f(t(i) + h/2, y(i) + (h/2)*k2);

k4 = f(t(i) + h, y(i) + h*k3);

y(i+1) = y(i) + h/6*(k1 + 2*k2 + 2*k3 + k4);

end

end

The function takes the function we want to evaluate, the bounds(a, b), the starting point, and the mesh size,

f = @(y,t) -t*y + 4*t/y;

a = 0;

b = 2.5;

h = .1;

y0 = 1;

[answer_rk, t_rk] = rungekutta(f, a, b, y0, h)

The result is a matrix filled with NaN values, this is because the first iteration evaluates k1 as zero. I can't tell if my error is a coding error or a misunderstanding of the Runge-Kutta implementation, could someone help?

Answer 1

[Torsten]

Change the order of the input arguments for f:

f = @(t,y) ...

instead of

f = @(y,t) ...

Answer 2

[Steven Lord]

You define f as:

f = @(y,t) -t*y + 4*t/y;

You call f as:

k1 = f(t(i), y(i));

Note the order of inputs in your definition and your call. Do you want f to be a function of y then t (as per the definition) or a function of t then y (as per the call)?