Runge-Kutta 4th order function error (Matrix dimensions must agree)

hello, i'm trying to make a Runge-Kutta 4th order function, and it works for some questions but i can't get it to work for this one, please check my code

close all; clc; clear all;

tspan = [0 10];

xi = [1;1];
params = [1 1]; %a = 1, b = 1 ; see function file

[t,x] = ode45(@hw3p4_func,tspan,xi,[],params);

y0 = x(1); y1 = x(2); ti = tspan(1) ; tf = tspan(2);


[t_rk, x_rk] = RK_4(@hw3p4_func, xi,y0,0.1);

plot(t,x(:,1))

%function to be solved
function dxdt = hw3p4_func(t,x, params)

y1 = x(1);
y2 = x(2);

a = params(1);
b = params(2);

dxdt = [a*y2;
        (-b*y1+(1-y1^2)*y2)];
end


%RK code
function [x,y] = RK_4(func, x,y0,delta_x)

y = zeros(1,length(x));
y(1) = y0;
n = length(x)-1;

for i = 1:n
     
    k1 = func(x(i),y(i));
    k2 = func(x(i)+.5.*delta_x,y(i)+.5.*k1.*delta_x);
    k3 = func(x(i)+.5.*delta_x,y(i)+.5.*k2.*delta_x);
    k4 = func(x(i)+    delta_x,y(i)+    k3.*delta_x);
    y(i+1) = y(i)+((k1+2.*k2+2.*k3+k4)./6).*delta_x;
    
end

Answer 1

[James Tursa]

Your RK_4 function is not set up to handle vector equations ... it is only set up to handle scalar equations. Also you are not passing params into RK_4, so the derivative function will not have that needed input.

So when you call it, you need to pass in an initial column vector and params in the function handle. E.g.,

[t_rk, x_rk] = RK_4(@(t,x)hw3p4_func(t,x,params), tspan, xi,0.1);

And the derivative function fixes would look something like this:

function [x,y] = RK_4(func, tspan,y0,delta_x)
n = round((tspan(2)-tspan(1))/delta_x)+1; % new calculation of n
x = linspace(tspan(1),tspan(2),n); % create a vector of x's from the endpoints with n elements
y = zeros(numel(y0),n); % changed to matrix of column vectors
y(:,1) = y0; % changed y to column vector
for i = 1:n-1 % changed top index on rhs to n-1
     
    k1 = func(x(i)            ,y(:,i)                ); % changed y to column vector
    k2 = func(x(i)+.5.*delta_x,y(:,i)+.5.*k1.*delta_x); % changed y to column vector
    k3 = func(x(i)+.5.*delta_x,y(:,i)+.5.*k2.*delta_x); % changed y to column vector
    k4 = func(x(i)+    delta_x,y(:,i)+    k3.*delta_x); % changed y to column vector
    y(:,i+1) = y(:,i)+((k1+2.*k2+2.*k3+k4)./6).*delta_x; % changed y to column vector
    
end
end

You could plot the results as (picking off 1st row since states are stored as column vectors):

plot(t_rk,x_rk(1,:));