function soln = sersol(ode, x, degree, initial, x0)

% Series solution of ordinary differential equation.
%
% SERSOL(ode, x, degree, initial, x0) returns the series
% solution up to the specified degree (at x = x0) of the
% ordinary differential equation ODE with independent
% symbolic variable x.  The fourth argument is a vector
% of initial conditions, starting with y(0), y'(0), ....
% For a nonsingular equation, the number of initial
% conditions must equal the order of the equation.  If
% the equation has a regular singular point at 0, then
% only one initial condition, representing the first
% coefficient in the series, should be given.  The fifth
% argument can be omitted, in which case x0 = 0 is used.
%
% The first argument, ode, is a symbolic expression or
% string containing a formal Taylor polynomial expansion
% of the differential equation with nonzero constant
% term.  For a nonsingular equation, ode should be the
% result of substituting the expression y = a0 +
% a1*(x - x0) + a2*(x - x0)^2 + ..., which can be created
% with the function FORMALSERIES, into the equation.  If
% the equation has a regular singular point at 0, then
% the Frobenius series y = (x - x0)^r*(a0 + a1*(x - x0)
% + a2*(x - x0)^2 + ...) (for appropriate r) should be
% substituted instead, and before using SERSOL, the
% result must be multiplied or divided by an appropriate
% power of (x - x0) and simplified to make it a
% polynomial that is nonzero at x0.  In this case, the
% output of SERSOL represents the analytic part of the
% series for y; it should be multiplied by (x - x0)^r to
% get the series solution.
%
% Example:
%
% syms x
% y = formalseries(x, 5);
% dy = diff(y, x);
% d2y = diff(dy, x);
% ode = collect(d2y - x*dy - y, x);
% sersol(ode, x, 5, [1, 0])

% Error checking.
if nargin < 4
  error('Not enough input arguments.')
end

% Initialization.
if nargin < 5
  x0 = 0;
end
ode = sym(ode);
x = sym(x);
x0 = sym(x0);
len = length(initial);

% Set coefficients determined by the initial conditions.
for k = 0:(len-1)
  vars{k+1} = ['a' num2str(k)];
  vals{k+1} = sym(initial(k+1));
end

% Create the system of equations.
for k = len:degree
  eqn = subs(diff(ode, x, k - len), x, x0);
  eqn = subs(eqn, {vars{1:len}}, {vals{1:len}});
  eqns{k+1} = [char(eqn) '=0'];
  vars{k+1} = ['a' num2str(k)];
end

% Solve the equations for the unknown coefficients.
range = (len+1):(degree+1);
sol = solve(eqns{range}, vars{range});
for k = range
  vals{k} = sol.(vars{k});
end

% Substitute the coefficients into a formal series.
soln = subs(formalseries(x, degree, x0), vars, vals);