t = linspace(0,20,1000);
terms = 100;
a0 = 1;
a1 = 0;

a(1) = 0;
a(2) = 0;
a(3) = -1/(3*2)*a0
for n=2:terms
    a(n+2) = -1/((n+2)*(n+1))*a(n-1)
end

soln = ones(size(t),1)*a0;
for n=1:terms
   soln = soln + a(n)*t.^n 
end

plot(t,soln);
