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

a0 = x0;
soln = a0*t.^0;
a(1) = xd0;
soln = soln + a(1)*t.^1;
a(2) = 0;
soln = soln + a(2)*t.^2;
n = 0;
a(3) = -1/(2*(n+3)*(n+2))*a0;
soln = soln + a(3)*t.^3;
for n=1:100
   a(n+3) = -1/(2*(n+3)*(n+2))*a(n); 
   soln = soln + a(n+3)*t.^(n+3);
end
plot(t,soln,'m');