% integrates diff eqq for 313 EA sequence, given omega and initial
% conditions for the angles
close all
clc
clear all

omega = [0.01; 0; 0.01];%rad/s angular velocity
Y0 = [pi/3; pi/2; pi/4];%rad initial EA

tf = 100; %s 
TSPAN = [0 tf]; %duration of integration in seconds
OPTIONS = odeset('Maxstep', 1);
[TOUT,YOUT] = ode45(@(t,u) EA313diffEq(t, u, 10*sin(t)*omega),TSPAN,Y0,OPTIONS);
plot(TOUT,YOUT*180/pi)
grid on
hold on
xlabel('time [s]'); ylabel('EA [deg]')
legend('\Omega', 'i', '\omega')