Two things:
1 - By defining two different functions for y1 and y2 you decouple the system, which is the source of the main problem.
2 - After fixing point 1, I had to use a very small time step size to obtain an acceptable solution, indicating that the implementation of an adaptive step size might be beneficial to the efficiency of the code.
%% Initialization
tspan = 100;
h = 0.0005;
fact = 1/h;
x = 0:h:tspan;
y = zeros(2,length(x));
y(:,1) = [0.2;0]; % Initial conditions
%% Runge Kutta
for i = 1:(length(x)-1) % calculation loop
k_1 = odeSystem(x(i),y(:,i));
k_2 = odeSystem(x(i)+0.5*h,y(:,i)-0.5.*h.*k_1);
k_3 = odeSystem((x(i)+0.5*h),(y(:,i)-0.5*h*k_2));
k_4 = odeSystem((x(i)+h),(y(:,i)-k_3*h));
y(:,i+1) = y(:,i)+(1/6)*(k_1 + 2*k_2 + 2*k_3 + k_4)*h; % main equation
end
figure
plot(x,y(1,:),'r',x,y(2,:),'b')
axis([0 100 -1 1])
title('Runge-Kutta')
function dydt = odeSystem(t,y)
dydt(1,1) = y(2);
dydt(2,1) = -y(1);
end
