h1=25; %mm
h2=12.5; %mm
b=12.5; %mm

L=300; %mm

E=210000; % N/mm^2

areaStart=h1*b; %Anfangsquerschnitt
areaEnd=h2*b; %Endquerschnitt

nElements=10000; %Elementzahl

F=20000;


areaX=@(x)areaStart-x/L*(areaStart-areaEnd);
lElement=L/nElements;

K_global_u=zeros(nElements+1);
tic

for ii_=1:1:nElements
    
    areaMean=areaX(lElement*(ii_-0.5));
    
    K_lokal=E*areaMean/lElement*[1 -1;
                                 -1 1];

    %lokale Steifigkeitsmatrix zu globaler hinzuaddieren
    K_global_u(ii_:ii_+1,ii_:ii_+1)=K_global_u(ii_:ii_+1,ii_:ii_+1)+K_lokal;
     
end


K_global=K_global_u(2:end,2:end);

%% Gleichungssystem naeherungsweise loesen 
q=zeros(nElements,1);
q(end)=F;

y=inv(K_global)*q;


%% Gleichungssystem exakt loesen
up=@(x,t)F/E*1/areaX(x);
[x,u]=ode45(up,[0 300],0);

toc


%% Loesung plotten
nPos=(1:1:nElements)'*lElement; %Knotenpositionen
% close all
figure
plot([0;nPos],[0;y],'bx--')
hold on
plot(x,u,'g-');

xlabel('Knotenposition [mm]')
ylabel('Verschiebung [mm]')

u(end)