```function [x,dx,ddx] = Newmark(dt,x0,v0,F,M,K,C,varargin)
% [x,dx,ddx] = Newmark(dt,x0,v0,F,M,K,C,varargin) solves numerically the
% equations of motion of a damped system
%
% INPUT
% F : vector  -- size: [1x N] -- Time series representinf the time history of the load.
% M : scalar  -- size: [1 x 1] -- Modal mass
% K : scalar  -- size: [1 x 1] -- Modal stifness
% C : scalar  -- size: [1 x 1] -- Modal damping
% dt : scalar  -- size: [1 x 1] -- time step
% x0 : scalar  -- size: [1 x 1] -- initial displacement
% v0 : scalar  -- size: [1 x 1] -- initial velocity
% Varargin:
%           alpha (by default equal to 1/4)
%           beta (by default equal to 1/2)
%
% OUTPUT
% x: time history of the system displacement response
% dx: time history of the system velocity response
% ddx: time history of the system acceleration response
%
% Example:
% t = linspace(0,100,1000);
% dt = median(diff(t));
% y5 = zeros(size(t));
% clear inputFun
% Y = [0,10]';
% M = 1;
% K = 1;
% C = 0.005;
% F = cos(t); % expression of the harmonic force
% [y2] = Newmark(dt,Y(1),Y(2),F,M,K,C,'alpha',1/4)
%
% author: E. Cheynet. University of Stavanger.  last updated: 31/12/2016

%%
% options: default values
p = inputParser();
p.CaseSensitive = false;
p.parse(varargin{:});

% shorthen the variables name
alpha = p.Results.alpha ;
beta = p.Results.beta;

% initial acceleration
a0 = M\(F(1)-C.*v0-K.*x0);

N = numel(F);
ddx  = [a0,zeros(1,N-1)];
dx  = [v0,zeros(1,N-1)];
x  = [x0,zeros(1,N-1)];

for ii=1:N-1,