% Differenzenverfahren, Beispiel auf den Seiten 260 bis 262,
% Koeffizientenmatrix wird als "Sparse Matrix" gespeichert

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 ;

nRBlinks  = 3 ; % 3 Werte fuer Beschreibung der Einspannung
nRBrechts = 7 ; % 4 Werte fuer Beschreibung des freien Randes
nspare = (n-4)*5+nRBlinks+nRBrechts ;    
AA     = zeros (nspare,1) ;
ipos   = zeros (nspare,1) ;     
jpos   = zeros (nspare,1) ; 
b      = zeros (n     ,1) ;
iMitte = round (nA * 0.5 + 3.1) ;

%Einspannung links:
ipos(1:3) = [1  2 2] ;
jpos(1:3) = [3  2 4] ; 
AA  (1:3) = [1 -1 1] ;

ii = nRBlinks+1 ;
for i=3:n-2                             % Matrix A:
    ipos(ii:ii+4) = i ;
    jpos(ii:ii+4) = i-2:i+2 ;
    if (i == iMitte)
      AA  (ii:ii+4) = [0 0 1 0 0] ;
    else
      AA  (ii:ii+4) = [1 -4 6 -4 1] ;
      b   (i)       = q0*h^4/EI     ;
    end
    ii=ii+5 ;
end

% Freier Rand rechts:
ipos(ii:ii+6) = [n-1 n-1 n-1 n-1  n   n   n ] ;
jpos(ii:ii+6) = [n-4 n-3 n-1  n  n-3 n-2 n-1] ;
AA  (ii:ii+6) = [ -1  2   -2  1   1   -2  1 ] ;

A = sparse (ipos,jpos,AA,n,n) ;

t1 = cputime ;
v = A \ b ;    % Loesen des Gleichungssystems (Berechnung der Durchbiegung v)
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')