% Differenzenverfahren, Beispiel auf den Seiten 260 bis 262,
% Ausnutzung der Bandstruktur der Koeffizientenmatrix

clear all

L  = 1 ;
q0 = 1 ;
EI = 1 ;
nA = 1000   ;       % Anzahl der Abschnitte (nur dieser Wert muss geaendert werden)
n  = nA + 5 ;       % Anzahl der Gleichungen
h  = L / nA ;

ibwr = 3 ;
ibwl = 4 ;
ibw  = ibwl + ibwr - 1 ;
nibw = n * ibw ;
iMitte = round (nA * 0.5 + 3.1) ;

A  = zeros (nibw,1) ;  % Nullvektor fuer Aufnahme der Bandmatrix
b  = zeros (n,1)    ;  % Nullvektor

A(1:2*ibw)           = [ 0  0  0  0  0  1 ...
                         0  0  0 -1  0  1 ] ;       % Randbedingungen links: Einspannung
A((n-2)*ibw+1:n*ibw) = [-1  2  0 -2  1  0 ...
                         1 -2  1  0  0  0 ] ;       % Randbedingungen rechts: Freier Rand

for i = 3:n-2                          
    A((i-1)*ibw+1:i*ibw) = [0 1 -4 6 -4 1] ;        % Standardgleichungen
    b(i)                 = q0*h^4/EI       ;
end

A((iMitte-1)*ibw+1:iMitte*ibw) = [0 0 0 1 0 0] ;        % Standardgleichungen
b(iMitte) = 0 ;

% Lösen des Gleichungssystems 
t1 = cputime ;
v = gabamp (A , b , n , ibwr , ibwl) ;              % Berechnung der Durchsenkung
ZeitGlSyst = cputime - t1

format long
vEnd    = v(n-2)
phiEnd  = (-v(nA+2) + v(nA+4)) / (2*h)
Mblinks = -EI/h^2 * (v(2) - 2*v(3) + v(4))

for i = 3:n-2
    Mb(i) = - EI/h^2 * (v(i-1) - 2*v(i) + v(i+1)) ; 
end

z = 0 : h : L ;
subplot (2 , 1 , 1) ; plot (z , v (3:n-2)) , axis ij , title ('Verschiebung')
subplot (2 , 1 , 2) ; plot (z , Mb(3:n-2)) , title ('Biegemoment')