% ode21.m   second oder ode with initial conditions
function ode21
  % SOLVE  d^2x/d^2t + 5*dx/dt - 4 x = sin(10 t)
  % initial conditions: x(0) = 0, x'(0)=0
  t=0:0.001:3;   % time scale

  initial_x    = 0;
  initial_dxdt = 0;
  [t,x]=ode45( @rhs, t, [initial_x initial_dxdt] );

  plot(t,x(:,1),t,x(:,2));
  title('solve d^2x/d^2t + 5*dx/dt - 4 x = sin(10 t)');
  xlabel('t');
  ylabel('x');
  legend('dx_1','dx_2');

  function dxdt=rhs(t,x)
    dxdt_1 = x(2);
    dxdt_2 = -5*x(2) + 4*x(1) + sin(10*t);
    dxdt=[dxdt_1; dxdt_2];
  end % rhs
end % ode21
