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

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

for n=2:100
   a(n+2) = -1/(3*(n+1)*(n+2))*(n*(n-1)+2)*a(n); 
   soln = soln + a(n+2)*t.^(n+2);
end
plot(t,soln,'g');
axis([0 5 -2 8])