这是一个常规问题, 给出代码参考吧:
m=1; b=1; r=1; g=9.8; w=1; F0=1; % 设为自己需要的数.
fun=@(t,y) [y(2);b/m*y(2)-g/r*sin(y(1))+F0*cos(y(3))/m/r;w];
[t,y]=ode45(fun,[0,10],[1,1,1]);
plot3(y(:,1),y(:,2),y(:,3))
grid on
xlabel \theta
ylabel y
zlabel z