function a=DEQseries(t,y,n);
% DEQSERIES taylor series for a particular differential equation
%   a=DEQseries(t,y,n); computes n coefficients of the Taylor series
%   y(t+dt)=sum_{k=0} a_{k+1} dt^k of the solution of the differential
%   equation y'=t^2+y^2 using the power series ansatz

a=zeros(1,n);
a(1)=y;
a(2)=t^2+a(1)^2;
a(3)=t+a(1)*a(2);
a(4)=(1+2*a(1)*a(3)+a(2)^2)/3;
for k=4:n-1
   a(k+1)=a(1:k)*a(k:-1:1)'/(k);
end