function f=quicksort(x,left,right)
if left<right[i,x]=Division(x,left,right);x=quicksort(x,left,i-1);x=quicksort(x,i+1,right);
end
f=x;

function [i,x]=Division(x,left,right)
base=x(left,2);
while left<rightwhile left<right&x(right,2)>=baseright=right-1;endc=x(left,2);d=x(right,2);c1=x(left,1);d1=x(right,1);c2=x(left,3);d2=x(right,3);c3=x(left,4);d3=x(right,4);c5=x(left,5);d5=x(right,5);x(left,2)=d;x(right,2)=c;x(left,1)=d1;x(right,1)=c1;x(left,3)=d2;x(right,3)=c2;x(left,4)=d3;x(right,4)=c3;x(left,5)=d5;x(right,5)=c5;
%     x(left,1)=x(right,1);while left<right&x(left,2)<=baseleft=left+1;endc=x(left,2);d=x(right,2);c1=x(left,1);d1=x(right,1);c2=x(left,3);d2=x(right,3);c3=x(left,4);d3=x(right,4);c5=x(left,5);d5=x(right,5);x(left,2)=d;x(right,2)=c;x(left,1)=d1;x(right,1)=c1;x(left,3)=d2;x(right,3)=c2;x(left,4)=d3;x(right,4)=c3;x(left,5)=d5;x(right,5)=c5;
%     x(right,1)=x(left,1);
end
i=left;

function [xstar,index,it]=bisect(fun,a,b,ep)
if nargin<4 ep=1e-5;end
fa=feval(fun,a);
fb=feval(fun,b);
if fa*fb>0xstar=[fa,fb];index=0;it=0;return
end
k=0;
while abs(b-a)/2>epx=(a+b)/2;fx=feval(fun,x);if fx*fa<0b=x;fb=fx;elsea=x;fa=fx;endk=k+1;
end
xstar=(a+b)/2;index=1;it=k;

function chi2plot(X)
dd=[];
p=[];
[M,N]=size(X);
MEAN=mean(X);
SS_1=inv(cov(X));
for byk=1:M;DD=(X(byk,:)-MEAN)*SS_1*(X(byk,:)-MEAN)';dd=[dd,DD];pp=(byk-0.5)/M;p=[p,pp];
end
dd=sort(dd)'
xx=chi2inv(p,N)'
plot(xx,dd,'+'),lsline
xlabel('chi2quantitle')
ylabel('Sample generalized diatance')
title('chi2plot')

function [x,y]=Euler_correct(fun,a,b,n,y0)
%改进的Euler公式，其中
%fun为一阶微分方程的函数
%a,b为求解区间的左右端点
%n为等分区间；
%y0为初始条件
x=zeros(1,n+1);y=zeros(1,n+1);
h=(b-a)/n;
x(1)=a;y(1)=y0;
for k=1:nx(k+1)=x(k)+h;y0=y(k)+h*feval(fun,x(k),y(k));y(k+1)=y(k)+h/2*(feval(fun,x(k),y(k))+feval(fun,x(k+1),y0));
end