There is no 'f' term in this system of ODE. Following shows how to solve this using ode45
t = linspace(0, 10, 1000); % solve the equation for t in [0, 10]
ic = [0; 0; 0; 0]; % initial condition
[t, y] = ode45(@(t, x) odeFun(t, x), t, ic);
subplot(2,2,1)
plot(t, y(:,1))
title('x1')
subplot(2,2,2)
plot(t, y(:,2))
title('x2')
subplot(2,2,3)
plot(t, y(:,3))
title('x1dot')
subplot(2,2,4)
plot(t, y(:,4))
title('x2dot')
function dxdt = odeFun(t, x)
% x(1)=>x1, x(2)=>x2, x(3)=>x1', x(4)=>x2'
dxdt = zeros(4, 1);
dxdt(1) = x(3);
dxdt(2) = x(4);
dxdt(3) = 1/23.3*(5926*x(3) - 4566*x(4) + 4*x(1) - 0.39*x(2));
dxdt(3) = 1/240*(-4566*x(3) + 11133*x(4) - 0.39*x(1) + 1.34*x(2) - 1);
end
