[y,t,erm]=dedit2(y,t,erm);

function [y,t,erm]=dedit2(y,t,erm);
persistent a1 a2 alph1 alph2 er ex firstCall i j k ng yt ; if isempty(firstCall),firstCall=1;end; 

y_shape=size(y);y=reshape(y,1,[]);
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(ng), ng=0; end;
if isempty(alph1), alph1=0; end;
if isempty(alph2), alph2=0; end;
if isempty(a1), a1=0; end;
if isempty(a2), a2=0; end;
if isempty(er), er=0; end;
if isempty(ex), ex=0; end;
if isempty(yt), yt=0; end;
if firstCall,   alph1=[1.0d0];  end;
if firstCall,  alph2=[1.0d0];  end;
if firstCall,  ng=[5];  end;
firstCall=0;
erm = 0.0d0;
if( t==0.0d0 )
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
ex = 0.0d0;
if( t<=30.0d0 )
ex = exp(-2.0d0.*t);
end;
a2 = 1.0d0;
for j = 1 : ng;
a1 = 1.0d0;
for i = 1 : ng;
k = fix(i +(j-1).*ng);
yt = t.^(i+j-2).*ex.*a1.*a2;
er = abs(y(k)-yt);
erm = max(erm,er);
a1 = a1.*alph1./i;
end; i = fix(ng+1);
a2 = a2.*alph2./j;
end; j = fix(ng+1);
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
end %subroutine dedit2