I have made some modifications to the shapes of the arrays in your code and changed the representation of the system of ODEs as a cell array of function handles to an array of function handles. The code below also demonstrates how "ode45" can be used to solve the problem. Both methods give the same result!
function testODE
f = @(t,y) [y(1)-4.*y(2);-y(1)+y(2)];
y0 = [1;0];
t0 = 0;
T = 1;
N = 11;
[Y,t] = RK42d(f, t0, T, y0, N);
figure;
hold on;
plot(t,Y(:,1),'r');
plot(t,Y(:,2),'b');
sol = ode45(f,[0,1],y0);
plot(sol.x,sol.y(1,:),'go');
plot(sol.x,sol.y(2,:),'ko');
legend('RK42d - y1','RK42d - y2','ode45 - y1','ode45 - y2');
end
function [Y, t] = RK42d(f, t0, T, y0, N)
h = (T - t0)/(N - 1);
Y = zeros(N,2);
t = linspace(t0,T,N);
Y(1,:) = y0;
for i = 1:(N-1)
y = Y(i,:)';
k1 = f(t(i),y);
k2 = f(t(i) +0.5*h, y +0.5*h*k1);
k3 = f(t(i) +0.5*h, y +0.5*h*k2);
k4 = f(t(i) +h , y +h*k3);
Y(i+1,:) = y + (h/6)*(k1+ 2.*k2 + 2*k3 + k4);
end
end
For more information on "ode45", please see
Best Answer