% Solve a first-order scalar ODE

syms x(t)
eqn1 = t*diff(x(t),t) == x(t) + 3*t^2*cos(2 * t)
%sol1 = dsolve(eqn1)
ic1 = x(pi) == -4
sol1 = simplify(dsolve(eqn1,ic1))

%
%% Solve a second-order ODE

syms y(t)
Dy(t) = diff(y(t),t); D2y(t) = diff(y(t),t,2);
eqn2 = D2y(t) + 4*Dy(t) + 13*y(t) == 18*exp(-2*t) + 136*cos(5*t)
sol2 = simplify(dsolve(eqn2))
%ic2 = [y(0) == 4, Dy(0) == -1];
%sol2(t) = simplify(dsolve(eqn2,ic2))
%fplot(sol2(t))

%
%% Solve a system of ODEs

syms x(t) y(t)
eqn3 = [diff(x(t),t) == y(t) + exp(-t), diff(y(t),t) == -2*x(t) - 3*y(t) + 4]
[sol3x,sol3y] = dsolve(eqn3);
x3 = simplify(sol3x)
y3 = simplify(sol3y)
%ic3 = [x(0) == 2, y(0) == -5]
%[sol3x,sol3y] = dsolve(eqn3,ic3);
%x3(t) = simplify(sol3x)
%y3(t) = simplify(sol3y)
%fplot(x3(t),'b-')
%hold on
%fplot(y3(t),'ro')
%hold off
%xlim([0 5])

%
%% Do partial fractions

syms s
f = (3*s-1)/(s^2+6*s+8)
fp = partfrac(f)

%
%% Solve ODEs with Heaviside and Dirac delta functions

syms y(t)
sympref('HeavisideAtOrigin',1);
Dy(t) = diff(y(t),t); D2y(t) = diff(y(t),t,2);
eqn4 = D2y(t) + 7*Dy(t) + 10*y(t) == 3*heaviside(t-1) - 5*dirac(t-2)
%sol4 = simplify(dsolve(eqn4,'IgnoreAnalyticConstraints',false))
ic4 = [y(0) == 0, Dy(0) == 0];
sol4(t) = simplify(dsolve(eqn4,ic4,'IgnoreAnalyticConstraints',false))
fplot(sol4(t),[0 6])