[t,y]=ode45(@project,[0 10],[0;0;0;0]) ;
plot(t,y(:,1))
figure
plot(t,y(:,2),'r')
figure
plot(t,y(:,3),'g')
figure
plot(t,y(:,4),'m')
function dydt= project(t,y)
M=5; m=0.5; R=0.5; tau=10*sin(100*t); F=10*cos(10*t) ; ddx=2;
dydt=zeros(4,1);
dydt(1)=y(3);
dydt(2)=y(4);
dydt(3)=-(M*R^2*y(3) - 2*tau + 2*R^2*y(3)*m + 2*R*ddx*sin(y(1)))/(2*R*cos(y(1)));
dydt(4)=(M^2*R^2*y(3) + 3*M*R^2*y(3)*m + 2*ddx*sin(y(1))*M*R - 2*tau*M + 2*cos(y(1))*sin(y(1))*R^2*y(3)^2 + 2*R^2*y(3)*m^2 + 2*ddx*sin(y(1))*R*m + 2*F*cos(y(1))*R - 2*tau*m)/(2*R^2*cos(y(1))^2);
end
Best Answer