%% Problems from the Series Solutions Chapter without the use of MupAD or dsolve
% Suppose we want to solve the differential equation \mathrm{y''}\left(x\right)-x\;\mathrm{y'}\left(x\right)-y\left(x\right),\mathrm{y'}\left(0\right)=1,y\left(0\right)=2
% by series methods, say keeping terms up to degree 5.  We can proceed as follows.  From the initial conditions, the constant term is 2 and the coefficiet of x is 1.  We let a2 a3 a4 a5 represent the other coefficients.

syms x a2 a3 a4 a5
y = sum([2 1 a2 a3 a4 a5].*x.^(sym([0:5])))

%
% The next step is to solve for the unknown coefficients using the differential equation.

LHS = coeffs(diff(y, x, 2) - x * diff(y, x) - y, x)

%

[c2 c3 c4 c5] = solve(LHS(1:4) == sym(zeros(1, 4)), [a2 a3 a4 a5])

%
% Thus our solution is

subs(y, [a2 a3 a4 a5], [c2 c3 c4 c5])

%
% Now we'll solve a nonlinear first order equation y'(x) = y(x)^2 + x, y(0) = -3, using a vector of symbolic variables for the coefficients.

A = sym('A', [1, 10]);
f = -3 + A*transpose(x.^[1:10])
C = coeffs(diff(f, x) - f^2 - x, x)

%
% Let's say we want a series solution to order 6.

B = solve(C(1:6), A(1:6))

%
% So our approximate solution is:

-3 + [B.A1, B.A2, B.A3, B.A4, B.A5, B.A6]*transpose(x.^[1:6])

%
% Let's check this against the solution of the differential equation computed by dsolve expanded out in series.

syms y(x)
fsol = simplify(dsolve(diff(y, x) == y(x)^2 + x, y(0) == -3))
simplify(series(fsol, x))

% Same answer, basically!