solution:
omega = 1; alpha = 1;beta=1;
f=@(ts,theta)[omega + alpha*(sin((theta(2)-theta(1))));
omega + alpha*(sin(theta(1)-theta(2)))+beta*(sin(theta(3)-theta(2)));
omega + beta*(sin((theta(2)-theta(3))))];
range = [0,1];
thetainit = [0.1,0.3,0.2];
[range,theta] = ode45(f,range,thetainit);
h=plot(range,theta,'Linewidth', 2);
set(h,{'Color'},{'r'; 'g'; 'b'})
xlabel('Range');
ylabel('Theta');
legend('Theta1','Theta2','Theta3');
