function y=heatdrv(t,n,s,e)
% calcul de  y^(j)/(2j)! pour la fonction
%   y(t)=exp( -1/ (t/e*(1-t/e))^s ) qui est 1+1/s Gevrey-Romieu pour nu >0
% 
% entrees:
%    t(nt,1): un vecteur colonne de nt instants  entre 0 et e
%    n: un ordre de derivation
%    s: l'exposant 
%    e: nombre positif 
%  sorties: 
%    y(nt,n+1): une matrice 
%           y(i,j) est la valeur de  y^(j-1)/(2j-2)! en t(i)
%
% Pierre Rouchon
% mai 96 
% 
[nt,mt]=size(t);
if mt>1 
  error('derive: t doit etre un vecteur colonne')
end;

y=zeros(nt,n+1);
tau=t/e;
tau=max(0.001,tau);tau=min(0.999,tau); % evite pb num vers 0 et vers 1
p=tau.*(1-tau);
p1=(1-2*tau)/e;
p2=-2/e/e;

%  A  correspond a  ( t/e(1-t/e)^(s+1) ) ^(j) / (j)!
A=zeros(nt,n+1);
A(:,1)=p.^(s+1);
A(:,2)=(s+1)*p1.*A(:,1)./p;
for i=2:n;
  A(:,i+1)= ( (s+2-i)*p1.*A(:,i) + (s+2-i/2)*p2*A(:,i-1) ) ./  (i*p);  
end
% 
y(:,1)=exp( -1./(p.^s) );
y(:,2)=s*p1.*y(:,1)./A(:,1)/2;
for i=2:n;
   f=s*( p1.*y(:,i)+p2*y(:,i-1)/(2*i-3)/2 ) / (2*i)/(2*i-1) ;
   x=-1;
   for k=1:i-1
     x=x*(i-k)/(2*i-2*k+2)/(2*i-2*k+1);
     f=f+x*A(:,k+1).*y(:,i-k+1);
   end;
   y(:,i+1)=f./A(:,1);
end