[sunjq]

请试试
function dx=myf(t,x)
dx=zeros(3,1)
dx(1)=1;
dx(2)=1/(x(3)-x(1));
dx(3)=-1/x(2)+1;

[t,y]=ode45('myf',[0 1],[0 x0 y0]);