close all
clear all

global k;

k = 1;

options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-4 1e-4]);
[T,Y] = ode45(@rhs,[0 30],[0 0 0 0],options);

figure
plot(T,Y(:,1),'-',T,Y(:,2),'-.')
figure
plot(T,Y(:,3),'-',T,Y(:,4),'-.')

F21 = k*(Y(:,3)-Y(:,4));
F12 = -k*(Y(:,3)-Y(:,4));
figure
plot(T,F21,'-',T,F12,'-.')