[t,y] = ode45(@series2rhs,[0 5],[1 2]);
plot(t,y(:,1));
hold on;

a0 = 1;
a(1) = 2;
a(2) = -a0/2;
soln = a0*t.^0+ a(1)*t + a(2)*t.^2;

for n=1:50
   a(n+2) = -1/((n+1)*(n+2))*(n*(n+1)*a(n+1) + a(n)); 
   soln = soln + a(n+2)*t.^(n+2);
end
plot(t,soln,'k');
axis([0 5 -2 4])