clc
close all
%% DATA
a = ; % left interval end
b = ; % right interval end
kmax = ; % maximum number of iterations 
edop = ; % admissible error
N = ; % number of FE
M = ; % number of load increments
%
%% INITIAL CALCULATIONS
% boundary conditions
ya = ; % exact solution at x=a (solution formula for x=a)
yb = ; % exact solution at x=b (solution formula for x=b)
% coefficients of initial solution (linear function y = A*x + B),
% satisfying conditions y(a) = ya and y(b) = yb
A = ; % coefficient A (formula 3.8)
B = ; % coefficient B (formula 3.8)
[A B]
% plots of the exact and initial solutions
figure(1);
x =  % "dense" vector of points from a to b
plot(,,'k-'); % plot of the exact solution
hold on
grid on
plot(,,'r:'); % plot of the initial solution
%
%% FE CALCULATIONS
syms x
f = ; % right-hand side function
n = ; % number of nodes
h = ; % length of FE (interval length by number of FE)
X = ; % nodes coordinates (sequence from a by h to b)

Y_pop = ; % values of initial solution at nodes (formula 3.6 for x = X)
Y_pop = Y_pop';
for j =  % loop over load increments(from 1 to M)
    disp(['Load increment no = ',num2str(j),'/',num2str(M)])
    for k =  % loop over N-R iterations (from 1 to kmax)
        KT = ; % fill global tangent stiffness matrix with zeros [n x n] 
        R = ; % fill global residual vector with zeros [n x 1]
        for i =  % loop over FE (from 1 to N)
            x1 = ; % coordinate of the first node of i-th FE
            x2 = ; % coordinate of the second node of i-th FE
            N1 = ; % first linear shape function
            N2 = ; % druga linear shape function
            Ye = ; % extract "i" and "i+1" elements from Y_pop
            Ne = ; % vector of shape functions [1x2]
            Be = ; % vector of shape functions' derivatives [1x2]
            KTe = ; % local tangent stiffness matrix (formula 3.2)
            KTe = eval(KTe); 
            Re = ; % local residual vector (formula 2.6)
            Re = eval(Re);
            edof = ; % vector [1x3] with element number and its d.o.f.
            [KT,R] = assem(  ; % assembly of stiffness matrix and residual vector
        end
        bc = ; % [2 x 2] vector of fixed d.o.f.
        % (first row: number of first node and 0 value; second row: number of last nodes and zero value]
        dY = solveq( ; % solution of SLAE: KT*dY = -R (determination of solution increment)        
        Y = ; % current solution (sum of previous solution and solution increment)        
        e = ; % error norm
        disp(['iteration no = ',num2str(k),'/',num2str(kmax),' e = ',num2str(e)]);
        Y_pop = Y;
        if  % brake-off test
            break;
        end
    end
end

% plot of FE solution
plot(X,Y,'ro-');
legend('exact solution','initial solution','FE solution');