y=0:10:750000;
u0=9.12;
v=(1/u0).*y;
x=v.^.2;
p=gammainc(x,5);
i=-.5.*p+.5;
subplot(222);
%plot(y,i)

t=0:.05:10;
q=.5*erfc(t/(2^.5));
tt=2^.5*erfinv(1-2*i);
plot(tt,y)
z=0:-10:-750000;vv=(-1/u0).*z;
xx=vv.^.2;p1=gammainc(xx,5);ff=-2^.5*erfinv(p1);
subplot(223)
plot(ff,z)
a11=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1];[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0];n=21
a22=zeros(n-2, n);
a=[a11;a22];
for i=3:n
   for j=1:n-1
      a(i,j)=a(i-1,j+1)-(i-2)*a(i-2,j);
      a(i,n)=-(i-2)*a(i-2,n);
   end
   end
va=exp(-.5.*(tt.^2));
for i=1:n
pp=a(i,:);
f=polyval(pp,tt.*2.^(-.5));
d=va.*f.*y;

s1(i)=0;

for k=2:length(tt)
s1(i)=s1(i)+.5*(tt(k)-tt(k-1))*(d(k)+d(k-1));
end
end
va1=exp(-.5.*(ff.^2));
for s=1:n
pp1=a(s,:);
f1=polyval(pp1,ff.*2.^(-.5));
d1=va1.*f1.*z;
s2(s)=0;
for q=2:length(ff)
s2(s)=s2(s)+.5*(ff(q-1)-ff(q))*(d1(q)+d1(q-1));
end
end
s=s1+s2;

