Unable to perform assignment because the size of the left side is 1-by-2 and the size of the right side is 2-by-2. (Line 32)

function [ t,y ] = rk4sys( f,tspan,y0,h )
%function [ t,y ] = rk4sys( f,tspan,y0,h )
%Uses RK4 method to solve the system ode y1'=p(t,y1,y2) y1(t0)=y1_0
%y2'=q(t,y1,y2) y2(t0)=y2_0
%on the interval tspan=[t0, tend] using time-step h.
%f=[p;q] and is formatted as required by ode45
%y0=[y1_0, y2_0]
%y is the output vector approximating [y1(t), y2(t)] at times stored in t.


%calculate n the number of time-steps needed to reach tend.
n=ceil((tspan(2)-tspan(1))/h);

%Store initial conditions in t and y.
y(1,1)=y0(1);
y(1,2)=y0(2);
t(1)=tspan(1);

%RK4 Iteration Loop
for i=1:n
    
    %If on last time-step, redefine h to make final t value = tend. 
    if i==n
        h=tspan(2)-t(i); 
    end
    
    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*((k1+2*k2+2*k3+k4)/6);%RK4 iteration
    t(i+1,1)=t(i)+h;
end

figure(1)
plot(t,y(:,1),t,y(:,2))
legend('y1','y2')

figure(2)
plot(y(:,1),y(:,2))
xlabel('y1');
ylabel('y2');

[Riley McFerran]

My issue is on line 32, or: y(i+1,:)=y(i,:)+h*((k1+2*k2+2*k3+k4)/6);%RK4 iteration

I'm not sure how to fix it

[Riley McFerran]

This is the code I typed into the command window. After entering the second line, is when I get the error

pp=@(t,x)[.23*x(1)-.0133*x(1)*x(2); -.4*x(2)+.0004*x(1)*x(2)]

[t,x]=rk4sys(pp,[0,65],[550,20],0.3)

Answer 1

[VBBV]

pp=@(t,x)[.23*x(1)-.0133*x(1)*x(2); -.4*x(2)+.0004*x(1)*x(2)]

pp = function_handle with value:
    @(t,x)[.23*x(1)-.0133*x(1)*x(2);-.4*x(2)+.0004*x(1)*x(2)]

[t,y]=rk4sys(pp,[0,65],[550,20],0.3)

n = 217

t = 1×218

0 0.3000 0.6000 0.9000 1.2000 1.5000 1.8000 2.1000 2.4000 2.7000 3.0000 3.3000 3.6000 3.9000 4.2000 4.5000 4.8000 5.1000 5.4000 5.7000 6.0000 6.3000 6.6000 6.9000 7.2000 7.5000 7.8000 8.1000 8.4000 8.7000

y = 218×2

550.0000 20.0000 545.2491 18.9428 542.7727 17.9337 542.4326 16.9756 544.1123 16.0699 547.7151 15.2174 553.1617 14.4179 560.3886 13.6709 569.3460 12.9751 579.9967 12.3293

figure(1)
plot(t,y(:,1),t,y(:,2))
legend('y1','y2')

figure(2)
plot(y(:,1),y(:,2))
xlabel('y1');
ylabel('y2');

function [ t,y ] = rk4sys( f,tspan,y0,h )
%function [ t,y ] = rk4sys( f,tspan,y0,h )
%Uses RK4 method to solve the system ode y1'=p(t,y1,y2) y1(t0)=y1_0
%y2'=q(t,y1,y2) y2(t0)=y2_0
%on the interval tspan=[t0, tend] using time-step h.
%f=[p;q] and is formatted as required by ode45
%y0=[y1_0, y2_0]
%y is the output vector approximating [y1(t), y2(t)] at times stored in t.


%calculate n the number of time-steps needed to reach tend.
n=ceil((tspan(2)-tspan(1))/h)

%Store initial conditions in t and y.
y(1,1)=y0(1);
y(1,2)=y0(2);
t(1)=tspan(1);

%RK4 Iteration Loop
for i=1:n
    
    %If on last time-step, redefine h to make final t value = tend. 
    if i==n
        h=tspan(2)-t(i); 
    end
    
    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).';  % transpose this 
    k4=f(t(i)+h,y(i,:)+h*k3).'; % transpose this too 
    
    y(i+1,:)=y(i,:)+h*((k1+2*k2+2*k3+k4)/6);%RK4 iteration
    t(i+1)=t(i)+h;
end
% t
% y
end 

You need to transpose k3 and k4 too in function rk4sys

[Riley McFerran]

Thank you so much! Just overlooked that

[VBBV]

Welcome, please accept the answer if it resolved your problem, :)