%clear all;
%close all;
x0 = 1;
xd0 = 2;
[t,y] = ode45(@varcoef3rhs,[0 10],[x0 xd0]);
plot(t,y(:,1));
hold on;

a0 = x0;
soln = a0;
a(1) = xd0;
soln = soln + a(1)*t;
a(2) = -10/18*a0;
soln = soln + a(2)*t.^2;
a(3) = -10/54*a(1);
soln = soln + a(3)*t.^3;
for n=2:100
   a(n+2) = -((n*(n-1)+10)/(9*(n+2)*(n+1)))*a(n); 
   soln = soln + a(n+2)*t.^(n+2);
end
plot(t,soln,'k');
axis([0 10 -3 3])