clear

# Given the 2nd order DE x" = f(t,x,x')
# Let  x1 = x
# Let  x2 = x'
# Then x1' = x2
# And  x2' = f(t,x,x')
# so we have converted 1 2nd order DE into 2 1st order DEs
# and we can numerically integrate them 
# with Octave's built-in lsode function

function accel = f(t,x,v)
accel = (1-v) + sin(3*t);
endfunction

function xdot = g(x,t)
xdot(1) = x(2);
xdot(2) = f(t,x(1),x(2));
endfunction

xo = 0; # initial position
vo = 0; # initial velocity

# Using lsode, with initial conditions x(1)=initial_position 
# and x(2)=initial_velocity on t=[0,15] with 100 points:

x = lsode ("g", [xo;vo], (t = linspace (0, 2, 50)'));

figure(1); #position
plot(t,x(:,1))

figure(2); #velocity
plot(t,x(:,2))