Hey all, I've been trying to debug a model I've created of a 4DOF geared system within Simulink by reproducing the results in MATLAB going from a 2DOF upwards to verify the dynamics are acting correctly. At the moment I am purely trying to model the stiffness between two gears, my equations are simply: From my reduced Simulink model, applying an initial displacement of the pinion of 1e-3m I get this response (sorry its hard to make out): Obviously from the figure above, you can see the equal amplitude oscillations with opposite direction (as these gears are not on the same shaft). I zoomed into the plot and verified that the initial condition of the gear is what I set it to be. However, I've written the following MATLAB code to try and recreate the above, but to no avail. I'm not sure where I am going wrong here.. Mainscript Code (Text): clear all; clc; %% Defining global variables global r_p r_g I_p I_g k_mb % Assignment of ODE variable values: r_p = 0.1297; % Pinion radius (m) r_g = 0.1297; % Gear radius (m) I_p = 0.0025; % Pinion inertia (kg.m^2) I_g = 0.0025; % Gear inertia (kg.m^2) k_mb = 100; %========================================================================== %% initial coditions for ode solver %========================================================================== z0 = [1e-3,0,0,0]; % Initial Conditions for ODE45. t_span = [0 1]; % Timespan 0 - 1 second opts = odeset('Reltol', 1e-3, 'AbsTol', 1e-3); [t,z] = ode45(@ode2dof, t_span, z0, opts); %========================================================================== %% Post Processing %========================================================================== % Renaming the ode45 outputs for easier usage. theta_2 = z(:,1); theta_2d = z(:,2); theta_3 = z(:,3); theta_3d = z(:,4); % Plotting subplot(2,1,1) plot(t,theta_2,t,theta_3), xlabel('Time (s)'), ... ylabel('Angular Displacement (rad)'),title('Angular Displacements'),... legend('theta 2','theta 3','Location', 'northwest'),grid; subplot(2,1,2) plot(t,theta_2d,t,theta_3d), xlabel('Time (s)'), ... ylabel('Angular Velocity (rad/s)'), title('Angular Velocities'),... legend('theta 2d', 'theta 3d','Location', 'northwest'),grid; ODE Script Code (Text): function zdot = ode2dof(t,z) % defining global variables global r_p r_g I_p I_g k_mb %% S.S.E: zdot(2) = ((k_mb*r_p)/I_p)*(z(3)*r_g-z(1)*r_p); % theta_2dd zdot(4) = ((k_mb*r_g)/I_g)*(z(1)*r_p-z(3)*r_g); % theta_3dd z(2) = zdot(1); z(4) = zdot(3); % Outputs zdot = zdot'; The output plot looks like this instead: Any help would be appreciated, as I'm lost at what's going on here even though its quite a simple system.