function sijierungekutta(X0,Y0,Z0,h,T)
dxdt=inline('-0.013*x-1000*x*y', 'x', 'y', 'z');

dydt=inline('-2500*y*z', 'x', 'y', 'z');

dzdt=inline('-0.013*x-1000*x*y-2500*y*z', 'x', 'y', 'z');

x0 = [1 1 0]';
h=0.0001;
T=0.0201;
t = [0:0.0001:0.0201]';
x1 = ones(T/h+1,3);
x1(1,:) = x0;
x2 = ones(T/h+1,3);
x2(1,:) = x0;
f1 = ones(T/h,3);
f2 = ones(T/h,3);
f3 = ones(T/h,3);
f4 = ones(T/h,3);
format long

for i=1:T/h

f1(i,1) = feval（dxdt,x1(i,1),x1(i,2),x1(i,3)); 

f1(i,2) = feval（dydt,x1(i,1),x1(i,2),x1(i,3)); 

f1(i,3) = feval（dzdt,x1(i,1),x1(i,2),x1(i,3));

f2(i,1) = feval（dxdt,x1(i,1)+h/2*f1(i,1),x1(i,2)+h/2*f1(i,2),x1(i,3)+h/2*f1(i,3));

f2(i,2) = feval（dydt,x1(i,1)+h/2*f1(i,1),x1(i,2)+h/2*f1(i,2),x1(i,3)+h/2*f1(i,3));

f2(i,3) = feval（dzdt,x1(i,1)+h/2*f1(i,1),x1(i,2)+h/2*f1(i,2),x1(i,3)+h/2*f1(i,3));

f3(i,1) = feval（dxdt,x1(i,1)+h/2*f2(i,1),x1(i,2)+h/2*f2(i,2),x1(i,3)+h/2*f2(i,3));

f3(i,2) = feval（dydt,x1(i,1)+h/2*f2(i,1),x1(i,2)+h/2*f2(i,2),x1(i,3)+h/2*f2(i,3));

f3(i,3) = feval（dzdt,x1(i,1)+h/2*f2(i,1),x1(i,2)+h/2*f2(i,2),x1(i,3)+h/2*f2(i,3));

f4(i,1) = feval（dxdt,x1(i,1)+h*f3(i,1),x1(i,2)+h*f3(i,2),x1(i,3)+h*f3(i,3));

f4(i,2) = feval（dydt,x1(i,1)+h*f3(i,1),x1(i,2)+h*f3(i,2),x1(i,3)+h*f3(i,3));

f4(i,3) = feval（dzdt,x1(i,1)+h*f3(i,1),x1(i,2)+h*f3(i,2),x1(i,3)+h*f3(i,3));

x2(i+1,:) = x1(i,:)+h/6*(f1(i,:)+2*f2(i,:)+2*f3(i,:)+f4(i,:));

    x1(i+1,:) = x2(i,:);
end
format long
fprintf('\n      x          y         z         t   \n');
for j=1:T/h+1
fprintf('  %7.5f  %8.7f  %8.7f  %8.7f  \n',t(j),x2(j,1),x2(j,2),x2(j,3));
    
end

subplot(1,3,1);
plot(x2(1:1:T/h+1,1),x2(1:1:T/h+1,2),'r');%x-y的二维投影
title('x-y');

subplot(1,3,2)
plot(x2(1:1:T/h+1,1),x2(1:1:T/h+1,3),'r');%x-z的二维投影
title('x-z');

subplot(1,3,3)
plot(x2(1:1:T/h+1,2),x2(1:1:T/h+1,3),'r');%y-z的二维投影
title('y-z')；

% plot3(x1(1:1:T/h+1,1),x1(1:1:T/h+1,2),x1(1:1:T/h+1,3),'r');%三维图像
end