Solving ODE with Matlab

%Parameters
m_d=5.0; R=0.5;I_d=m_d*(R^2)/2;
m=2.5;
f_n=1;k=m*(2*pi*f_n)^2;
a = 0.4;
%Initial Conditions
y(1)=0;theta(1)=0;
v_y(1) = 0;tau_0=8;
w(1)=0;
dt=.00001; t_final=10;
t=0:dt:t_final;
   
for i=1:length(t) 
    %First order equations
   dy(i)=v_y(i);
    dtheta(i)=w(i);
    dw(i)=(-m*a*y(i)*(w(i))^2-2*m*y(i)*w(i)*v_y(i)+k*a*y(i)+tau_0)/(I_d+m*(y(i))^2);
    dv_y(i)=-a*w(i)+y(i)*(w(i))^2-(k/m)*y(i);
    %Integrating the equations using Euler integration
    y(i+1)= y(i)+dy(i)*dt;
    theta(i+1)=theta(i)+dtheta(i)*dt;
    w(i+1)=w(i)+dw(i)*dt;
    v_y(i+1)=v_y(i)+dv_y(i)*dt;
end
 
 
figure(1);
plot(t,w(1:length(t)),'k');grid on
xlabel('Time, s');ylabel('Angular velocity, rad/s');
figure(2);
plot(t,y(1:length(t)),'k');grid on
xlabel('Times,s');ylabel('Position of mass,m');

Here is my code. The problem is tau is function of time if t<=0.5s tau=tau_0 else tau=0. How can I write that in this loop?.

[Jan]

Today I've formatted your code. Please use the buttons over the edit section to do this by your own in future questions.

Answer 1

[Jan]

t = 0:dt:t_final;
for i = 1:length(t)
    % if t<=0.5s tau=tau_0 else tau=0
    if t(i) <= 0.5
        tau = tau_0;
    else
        tau = 0;
    end
    ...
end

The conversion from English to Matlab is easy in this case.

[Anh Dang]

Thank you so much :D