function y = seriessolve(eqn, yvar, xvar, y0, order)
%SERIESSOLVE solves the 1st order differential equation eqn == 0 for yvar 
% in terms of xvar, with initial condition yvar = y0 at xvar = 0,
% out to the specified order.  syntax:
% sol = seriessolve(eqn, yvar, xvar, y0, order)
%
% example: seriessolve(diff(y) - y, y, t, 1, 3) yields
% t^3/6 + t^2/2 + t + 1, the beginning of the Taylor expansion of exp(t)
%
A = sym('A', [1, order]); B = sym('B', [1, order]);
% create vectors of symbolic coefficients
f = y0 + A*transpose(xvar.^(sym([1:order]))); % create Taylor polynomial
eqn1 = subs(eqn, yvar, f);
C = coeffs(eqn1, xvar); C = C(0); %setting t = 0 is to convert to constants
for j = 1:order
    B(j) = solve(C(j) == 0, A(j)); C = subs(C, A(j), B(j));
end
y = subs(f, A, B);
end
