Inconsistent results using ODE45

Hi Guys

I have 3 differential equations which I have solved using the ode45 solver.

The state representation and solving using ode45 are presented in the below code.

When assisnging initial conditions to only the X and Y component the results from the varioustime increments are the same.

Once I assign initial conditions to the "theta" component the results using various time increments are all different as presented by the image.

The results vary so much that if processes the signals to attain the frequency domain, the frequencies even differ.

Why is this occuring and how can I resolve the issue? I am not sure how to deal with the current issue.

clc;
clear;
r = 0.5;
L = 2;

k = 266*1000/2;
m = 3755;

SolverOptions = odeset('RelTol',1e-5,'AbsTol',1e-5);
T = 30;
[t1,y1] = ode45(@Statespace3dof2,[0:0.1:T],[0.1 0 r+0.1 0 deg2rad(5) 0],SolverOptions,r,L,k,m); 
[t2,y2] = ode45(@Statespace3dof2,[0:0.01:T],[0.1 0 r+0.1 0 deg2rad(5) 0],SolverOptions,r,L,k,m); 
[t3,y3] = ode45(@Statespace3dof2,[0:0.001:T],[0.1 0 r+0.1 0 deg2rad(5) 0],SolverOptions,r,L,k,m); 
[t4,y4] = ode45(@Statespace3dof2,[0:0.0001:T],[0.1 0 r+0.1 0 deg2rad(5) 0],SolverOptions,r,L,k,m);
[t5,y5] = ode45(@Statespace3dof2,[0:0.00001:T],[0.1 0 r+0.1 0 deg2rad(5) 0],SolverOptions,r,L,k,m); 

figure(1)
hold on
L1 = plot(t1,y1(:,1));
L2 = plot(t2,y2(:,1));
L3 = plot(t3,y3(:,1));
L4 = plot(t4,y4(:,1));
L5 = plot(t5,y5(:,1));
legend([L1 L2 L3 L4 L5],'0.1','0.01','0.001','0.0001','0.00001')
title('X component');
hold off

figure(2)
hold on
L1 = plot(t1,y1(:,3));
L2 = plot(t2,y2(:,3));
L3 = plot(t3,y3(:,3));
L4 = plot(t4,y4(:,3));
L5 = plot(t5,y5(:,3));
legend([L1 L2 L3 L4 L5],'0.1','0.01','0.001','0.0001','0.00001')
title('Y component');
hold off


figure(3)
hold on
L1 = plot(t1,y1(:,5));
L2 = plot(t2,y2(:,5));
L3 = plot(t3,y3(:,5));
L4 = plot(t4,y4(:,5));
L5 = plot(t5,y5(:,5));
legend([L1 L2 L3 L4 L5],'0.1','0.01','0.001','0.0001','0.00001')
title('theta component');
hold off

function xp = Statespace3dof2(t,X,r,L,k,m)

x = X(1);
x_dot = X(2);
y = X(3);
y_dot = X(4);
theta = X(5);
theta_dot = X(6);


A = L^2/4+L*x+x^2+y^2-0.5*L^2*cos(theta)-L*x*cos(theta)+0.25*L^2*cos(theta)^2-L*y*sin(theta)+0.25*L^2*sin(theta)^2;
B = L^2/4-L*x+x^2+y^2-0.5*L^2*cos(theta)+L*x*cos(theta)+0.25*L^2*cos(theta)^2+L*y*sin(theta)+0.25*L^2*sin(theta)^2;

%theta
one_line1 = 0.5*k*L^2*sin(theta) + 0.25*k*L^2*(-2*cos(theta)*sin(theta))+0.25*k*L^2*(2*sin(theta)*cos(theta));
one_line2 = -k*r*(0.5*L^2*sin(theta) + L*x*sin(theta) + 0.25*L^2*(-2*cos(theta)*sin(theta)) - L*y*cos(theta) + 0.25*L^2*(2*sin(theta)*cos(theta)))/(2*sqrt(A));
one_line3 = -k*r*(0.5*L^2*sin(theta) - L*x*sin(theta) + 0.25*L^2*(2*cos(theta)*sin(theta)) + L*y*cos(theta) + 0.25*L^2*(2*sin(theta)*cos(theta)))/(2*sqrt(B));
one = one_line1+one_line2+one_line3;
theta_dotdot = -4*one/(L^2*m);

%X
two = 2*k*x-k*r*(L+2*x-L*cos(theta))/(2*sqrt(A))-k*r*(-L+2*x+L*cos(theta))/(2*sqrt(B));
x_dotdot = -two/(m);

%Y
three = -9.81*m+2*k*y-k*r*(2*y-L*sin(theta))/(2*sqrt(A))-k*r*(2*y+L*sin(theta))/(2*sqrt(B));
y_dotdot = -three/(m);

xp = [X(2)
      x_dotdot;
      X(4);
      y_dotdot;
      X(6);
      theta_dotdot;];
end