y[t_]=x[t]/.DSolve[{x'[t]==t+x[t]^2,x[0]==1},x[t],t][[1]]
Plot[y[t],{t,0,20}]
t0=.1;
y[t0]
Do[Print[Normal[Series[y[t],{t,0,j}]]/.t->t0],{j,1,7}]

Clear[a,y,lhs,rhs,t]
y[t_]=1+Sum[a[k] t^k,{k,1,10}];
lhs[t_]=D[y[t],t];
rhs[t_]=Series[t+y[t]^2,{t,0,10}];
Do[a[k]=a[k]/.Solve[SeriesCoefficient[lhs[t]-rhs[t],{t,0,k-1}]==0,a[k]][[1]];
Print[a[k]]
,{k,1,9}]
y[0.1]

(* Euler's method *)
z=1.;t=0;
f[t_,z_]:=t+z^2;
h=10^(-3);
tf=0.1;
n=tf/h;
Do[z=z+h*f[t,z];t=t+h
,{j,1,n}]
z