How do you make a script which uses the third- order Runge- Kutta scheme to solve for y'

How do I make a function which uses

to solve

I think that the function should have the interface

function y = odesRK3(fun, t, y0)
% where y is a matrix conatining the solution,
% fun is the user-supplied function f(t,y),
% t is a row vector containing values of the independant variable t in
% ascending order,
% and y0 is a column vector of initial conditions


end

To test if the function works, im using the following code to call the function

f = @(t,y) [-2*y(1)+y(2); 2*y(1)-y(2)];
h = 0.2;
t = 0:h:1.4;
y0 = [9;3];
y = odesRK3(f, t, y0);
yExact = [4+5*exp(-3*t); 8-5*exp(-3*t)];
overallError = max(abs(y-yExact), [], 2)
plot(t, y, 'o--', t, yExact)
xlabel('t')
legend('RK3 y_1','RK3 y_2','Exact y_1','Exact y_2')

btw, im assuming that fun(t, y) returns a column vector given a scalar t and a column vector y

also, the solution at each point in the vector t must be stored in the columns of y, as in, y(: ,k) should be the solution at t(k)

Answer 1

[Torsten]

A class mate of yours seems to have the same problem:

https://de.mathworks.com/matlabcentral/answers/1827643-function-uses-the-third-order-runge-kutta-scheme?s_tid=srchtitle

Maybe you can exchange ideas.

And - together we are strong - here is the combined code:

f = @(t,y) [-2*y(1)+y(2); 2*y(1)-y(2)];
h = 0.2;
t = 0:h:1.4;
y0 = [9;3];
y = odesRK3(f, t, y0);
yExact = [4+5*exp(-3*t); 8-5*exp(-3*t)];
overallError = max(abs(y-yExact), [], 2)

overallError = 2×1

0.1760 0.1760

plot(t, y, 'o--', t, yExact)
xlabel('t')
legend('RK3 y_1','RK3 y_2','Exact y_1','Exact y_2')

function y = odesRK3(f, t, y0)
% where y is a matrix conatining the solution,
% fun is the user-supplied function f(t,y),
% t is a row vector containing values of the independant variable t in
% ascending order,
% and y0 is a column vector of initial conditions
n=length(t);
y=nan(length(y0),n);
y(:,1)=y0(:);
for k=1:n-1
    h=t(k+1)-t(k);
    F1=h*f(t(k),y(:,k));
    F2=h*f(t(k)+h/2,y(:,k)+(F1/2));
    F3=h*f(t(k)+0.75*h,y(:,k)+(0.75*F1));
    y(:,k+1)=y(:,k)+(2*F1+3*F2+4*F3)/9;
end
end