Use ode45() for estimating a numerical solution
IC = [0.1, 0.3, 0.2];
tspan = [0 1];
[t, thetas] = ode45(@odeFun, tspan, IC);
plot(t, thetas);
legend({'\theta_1', '\theta_2', '\theta_3'}, 'FontSize', 16, 'Location', 'best');
function dthetas = odeFun(t, thetas)
theta1 = thetas(1);
theta2 = thetas(2);
theta3 = thetas(3);
dtheta1 = 1+sin(theta2-theta1);
dtheta2 = 1+sin(theta1-theta2)+sin(theta3-theta2);
dtheta3 = 1+sin(theta2-theta3);
dthetas = [dtheta1; dtheta2; dtheta3];
end
