% Solve an IVP for the linear second order ODE

%    y'' + p(t) y' + q(t) y = g(t)
%    y(t_0)  = y0
%    y'(t_0) = y1

% Define coefficient functions
p = @(t) 0;
q = @(t) -1;
g = @(t) 0;

% Define initial conditions
y0 = 0.0;
y1 = 1.0;

% Define initial and final times
t0 = 0.0;
t1 = 2.0;

% Write as first order system x' = f(t,x)
% x = [x(1);x(2)] with x(1) = y, x(2) = y'
%    x'(1) = x(2)
%    x'(2) =  - q(t)*x(1) - p(t)*x(2) + g(t)

% Define right hand side in ODE
f = @(t,x) [x(2);-q(t)*x(1)-p(t)*x(2)+g(t)];

% Define initial value
xzero(1) = y0;
xzero(2) = y1;

% Solve ODE
[t x] = ode45(f,[t0,t1],xzero);

% Plot results
 plot(t,x(:,1),'b',t,x(:,2),'r'), xlabel('t'), ylabel('y, dy/dt')
% plot(t,x(:,1),'b'), xlabel('t'), ylabel('y')