clc clear tspan = [0:0.001:20]; y0 = [0; 0; 0; 0]; [t,y] = ode45(@(t,y) Q4d_forced_2DOF(t,y), tspan, y0); figure plot(t,y(:,1)) % y1=x1 xlabel('t') ylabel('x_1') figure plot(t,y(:,3)) % y3=x2 xlabel('t') ylabel('x_2')