nonlinear ode of a cart pole system in matlab

main.m
function dydt = cart_pole_ode(t, y, m, M, L, g, b, u)
% y = [x; dx; theta; dtheta]
% m = mass of the pole
% M = mass of the cart
% L = length of the pole
% g = gravitational acceleration
% b = friction coefficient
% u = control input (force)

x = y(1); dx = y(2); theta = y(3); dtheta = y(4);

s = sin(theta); c = cos(theta);
D = m*L*L*(M+m*(1-c^2)); % denominator for the ODEs

dxdt = dx;
dtheta_dt = dtheta;

dx_dt_dt = (1/D)*(-m^2*L^2*g*c*s + m*L^2*(m*L*dtheta^2*s-c*dx*dtheta) ...
              - m*L*c*b*dx + m*L*u);

dtheta_dt_dt = (1/D)*((M+m)*m*g*L*s - m*L*c*(m*L*dtheta^2*s - c*dx*dtheta) ...
              + m*L*c*b*dtheta - m*L*c*u);

dydt = [dx; dx_dt_dt; dtheta; dtheta_dt_dt];            
end
707 chars
26 lines

main.m
% System parameters
m = 0.1;    % mass of the pole
M = 1;      % mass of the cart
L = 0.5;    % length of the pole
g = 9.81;   % gravitational acceleration
b = 0.1;    % friction coefficient

% Initial conditions
x0 = 0;       % initial position of the cart
dx0 = 0;      % initial velocity of the cart
theta0 = pi;  % initial angle of the pole
dtheta0 = 0;  % initial angular velocity of the pole

y0 = [x0; dx0; theta0; dtheta0];

% Control input
u = 0;   % no control input

% Solve the ODEs
[t, y] = ode45(@(t,y) cart_pole_ode(t,y,m,M,L,g,b,u), [0 10], y0);
562 chars
21 lines

main.m
% Plot results
x = y(:,1); theta = y(:,3);

figure;
subplot(2,1,1);
plot(t, x);
xlabel('Time (s)'); ylabel('Position (m)');

subplot(2,1,2);
plot(t, theta);
xlabel('Time (s)'); ylabel('Angle (rad)');
200 chars
12 lines