y=0:50:4500000;pas=50
t0=clock;
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:-50:-4500000;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 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 0 0 0 0 0 0 0 0 0 1 0];n=30
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);
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);
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;
s11(1)=s1(1);
s22(1)=s2(1);
w=1;
for l=2:n
   w=1;
   for j=1:l-1
      w=w*j;
   end
   s11(l)=s1(l)/w;
   s22(l)=s2(l)/w;
end
ss=s11+s22;sf=ones(1,length(ss));
sf=4.*ss.*s1;
durata=etime(clock,t0)
for i=1:n
   sf1(i)=sf(n+1-i);
end
rg=-1:.05:1;
ry=polyval(sf1,rg);
ry1=(1/max(ry)).*ry;
plot(rg,ry1)
qq=max(tt)
dom=erf(qq)