clear all
clc
L=1;
I=1;
E=1;
QL=1;
ML=0;
q=@(x)0;
mu=1;
n=5;
a=0;
b=0;
bc=[1,1,a;1,2,b;n,3,ML;n,4,QL];
x=linspace(0,L,n);
hj=1;


[M,S]=systemmatrices(x,E,I,mu);

[Me,Se,r]=linearsystem(x,M,S,q,bc);



% Number of time steps.
nT = 1e3;
% Create figure object.
fig = figure;

% Create axes object.
ax = axes('Box', 'on', 'DataAspectRatio', [1, 1, 1]);

% Initialize empty plot of deflection and return handle.
handle_deflection = line('XData', [], 'YData', [], 'LineStyle', 'none', 'Marker','x', 'MarkerEdgeColor', 'r');

%...


% Solve.
y = Se\r;

dy=zeros(2*n+2,1);
ddy=zeros(2*n+2,1);

% Set x and y-coordinate limits based on x and y.
xlims = [1 10];
ylims = [1 10];
set(ax, 'XLim', xlims, 'YLim', ylims);

% Iteration for all time steps.
for t = 1:nT
% Check if figure exists, else break.
if ~ishandle(fig)
    break
end
% Get w, dw;
w = y(1:2:end-2);
dw =y(2:2:end-2);

% Update deflection plot.
set(handle_deflection, 'XData', x, 'YData', w)

% Do a single NEWMARK?iteration step.
[y, dy, dyy] = newmark(y, dy, ddy, Me, Se, r, hj);

end