I wish to simulate the output of a certain gear system I have. How the gear system looks isn't particularly important to the problem, I managed to get the differential equation needed from the mechanical system. Here is the code I have

% parameters
N2  = 90;
N1  = 36;
Jn1 = 0.5;
Jn2 = 0.8;
J2  = 2;
D   = 8;
K   = 5;
J   = (N2/N1)^2 * Jn1 + Jn2 + J2;

% define the system
sys = ss([0 1; -K/J -D/J], [0; N2/(N1*J)], [1 0], 0);

% initial state: (position, velocity) [rad; rad/s]
x0 = [0; 0];

% define the time span
t = linspace(0, 15, 10000)';

% define the input step
T1 = zeros(length(t), 1);
T1(t>=0) = 1;

% compute the system step response at once
theta1 = lsim(sys, T1, t, x0);

% compute the system response as aggregate of the forced and unforced
% temporal evolutions
theta2 = lsim(sys, T1, t, [0; 0]) + initial(sys, x0, t);

% plot results
figure('color', 'white');
hold on;
yyaxis left;
plot(t, T1, '-.', 'linewidth', 2);
ylabel('[N]');
yyaxis right;
plot(t, theta1, 'linewidth', 3);
plot(t, theta2, 'k--');
xlabel('t [s]');
ylabel('[rad]');
grid minor;
legend({'$T_1$', '$\theta_1$', '$\theta_2$'}, 'Interpreter', 'latex',...
       'location', 'southeast');
hold off;

This should work in generating a graph that shows the positions, my outputs, for a Heaviside/step input. My question is, how would I do this for a sine wave input. I figure I should have sin(w*t) instead of (t>=0), where w is my pulse frequency. Still, I can't seem to make this work. Any help would be really appreciated! :)

Test Run

These are the output results I currently get on MATLAB R2019b. As Luis' comment has suggested I have also declared a sinusoid as T1 to serve as the input. Currently not sure if this result is the expected output.

Code Snippet:

t = linspace(0, 15, 10000)';
f = 0.1; 
phi = 0; 
T1 = sin(2*pi*f*t + phi);

f → Frequency of sinusoidal input (0.1Hz in this example).
phi → Phase offset of sinusoidal input/initial phase (0 in this example).
t → Time vector dictating the samples of the sinusoid.
0 → Start time (0 seconds in this example).
15 → End time (15 seconds in this example).
10000 → Number of samples between the start time (0s) and end time (15s).

Implementation in Script:

% parameters
N2  = 90;
N1  = 36;
Jn1 = 0.5;
Jn2 = 0.8;
J2  = 2;
D   = 8;
K   = 5;
J   = (N2/N1)^2 * Jn1 + Jn2 + J2;

% define the system
sys = ss([0 1; -K/J -D/J], [0; N2/(N1*J)], [1 0], 0);

% initial state: (position, velocity) [rad; rad/s]
x0 = [0; 0];

% define the time span
t = linspace(0, 15, 10000)';

% define the input step
T1 = zeros(length(t), 1);
T1(t>=0) = 1;

f = 0.1; %Sinusoid frequency = 0.1Hz%
phi = 0; %Phase = 0%
T1 = sin(2*pi*f*t + phi);

% compute the system step response at once
theta1 = lsim(sys, T1, t, x0);

% compute the system response as aggregate of the forced and unforced
% temporal evolutions
theta2 = lsim(sys, T1, t, [0; 0]) + initial(sys, x0, t);

% plot results
figure('color', 'white');
hold on;
yyaxis left;
plot(t, T1, '-.', 'linewidth', 2);
ylabel('[N]');
yyaxis right;
plot(t, theta1, 'linewidth', 3);
plot(t, theta2, 'k--');
xlabel('t [s]');
ylabel('[rad]');
grid minor;
legend({'$T_1$', '$\theta_1$', '$\theta_2$'}, 'Interpreter', 'latex',...
       'location', 'southeast');
hold off;