function soln = sersol(ode, x, order, initial)
x = sym(x);
syms u;
u = formalseries(x,order);
du = diff(u);
d2u = diff(du);

syms y dy d2y;
ode = sym(ode);
ode = subs(ode,y,u);
ode = subs(ode,dy,du);
ode = subs(ode,d2y,d2u);
serode = collect(ode,x);

d = length(initial);
for i=0:d-1
    tmp = strcat('a', num2str(i));
    tmp2 = num2str(initial(i+1));
    eval(strcat(tmp, ' = ', tmp2, ';'));
end

for i=d:order
    eqn = subs(diff(serode, x, i-d), x, 0);
    eqns{i} = strcat(char(eqn), ' = 0');
    tmp = strcat('a', num2str(i));
    vars{i} = tmp;
end

sol = solve(eqns{d:order}, vars{d:order});
for i=d:order
    tmp = strcat('a', num2str(i));
    tmp2 = char(getfield(sol, tmp));
    eval(strcat(tmp, ' = ', tmp2, ';'));
end
soln = subs(formalseries(x, order));


