fun=@(t,x)[36*(x(2)-x(1));-x(1)*x(3)-16*x(1)+28*x(2)-x(4);... x(1)*x(2)-3*x(3);x(1)+0.2][t,X] = ode45(fun,[0:0.01:20],[0 0 0 0]);plot(t,X)
用ode函数，见帮助
可以用mathmatica直接算只要给出式子就能算
