% Example 2.3 Two nonlinear springs (Newton-Raphson method)
tol = 1.0e-5; iter = 0; c = 0;            % Tolerance and initialize parameters
u = [0; 0]; uold = u;                                        % Initial solution
f = [0; 100];                                            % Applied force vector
P = [300*u(1)^2+400*u(1)*u(2)-200*u(2)^2+150*u(1)-100*u(2) % Nonlinear function
     200*u(1)^2-400*u(1)*u(2)+200*u(2)^2-100*u(1)+100*u(2)];
R = f - P;                                                           % Residual
conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);                          % Convergence
fprintf('\n iter      u1      u2         conv       c');           % Print head
fprintf('\n %3d  %7.5f %7.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
while conv > tol && iter < 20               % Loop for Newton-Raphson iteration
 Kt = [600*u(1)+400*u(2)+150  400*(u(1)-u(2))-100             % Jacobian matrix
       400*(u(1)-u(2))-100    400*u(2)-400*u(1)+100];
 delu = Kt\R;                                    % Solve for solution increment
 u = uold + delu;                                             % Update solution
 P = [300*u(1)^2+400*u(1)*u(2)-200*u(2)^2+150*u(1)-100*u(2);  % Nonlinear func.
      200*u(1)^2-400*u(1)*u(2)+200*u(2)^2-100*u(1)+100*u(2)];
 R = f - P;                                                          % Residual
 conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);                         % Convergence
 c = abs(0.9-u(2))/abs(0.9-uold(2))^2;                   % Convergence constant
 uold = u;                                  % Store solution for next iteration
 iter = iter + 1;                                  % Increase iteration counter
 fprintf('\n %3d  %7.5f %7.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
end
%
% Example 2.4 Divergence of the Newton-Raphson method.
xdata=zeros(40,1); ydata=zeros(40,1);             % Store history of iterations
tol = 1.0e-5; iter = 0; c=0;              % Tolerance and initialize parameters
u = 0.5; uold = u;                                           % Initial solution
P = u+atan(5*u); R = -P;                        % Nonlinear function & residual
conv= R^2;                                                        % Convergence
xdata(1)=u; ydata(1)=P;                         % Initial solution and function
while conv > tol && iter < 20               % Loop for Newton-Raphson iteration
 Kt = 1+5*(cos(atan(5*u)))^2;                                 % Jacobian matrix
 delu = R/Kt;                                    % Solve for solution increment
 u = uold + delu;                                             % Update solution
 P = u+atan(5*u); R = -P;                       % Nonlinear function & residual
 conv= R^2;                                                       % Convergence
 uold = u;                              % Store solution for the next iteration
 iter = iter + 1;                                   % Increse iteration counter
 xdata(2*iter)=u;   ydata(2*iter)=0;             % Linearized solution location
 xdata(2*iter+1)=u; ydata(2*iter+1)=P;                    % Vertical error in P
end                                      % End of Newton-Raphson iteration loop
plot(xdata,ydata); hold on;                            % Plot iteration history
x=-1:0.1:1; y=x+atan(5*x); plot(x,y);                  % Plot the true function
%
% Example 2.5 Nonlinear springs (modified Newton-Raphson method)
tol = 1.0e-5; iter = 0; c = 0;            % Tolerance and initialize parameters
u = [0.3; 0.6]; uold = u;                                    % Initial solution
P = [300*u(1)^2+400*u(1)*u(2)-200*u(2)^2+150*u(1)-100*u(2) % Nonlinear function
     200*u(1)^2-400*u(1)*u(2)+200*u(2)^2-100*u(1)+100*u(2)];
f = [0; 100]; R = f - P;                           % Applied force and residual
conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);                          % Convergence
fprintf('\n iter      u1      u2         conv       c');           % Print head
fprintf('\n %3d  %7.5f %7.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
Kt = [600*u(1)+400*u(2)+150  -400*u(2)+400*u(1)-100           % Jacobian matrix
      400*u(1)-400*u(2)-100  400*u(2)-400*u(1)+100];
while conv > tol && iter < 20               % Loop for Newton-Raphson iteration
 delu = Kt\R;                                    % Solve for solution increment
 u = uold + delu;                                             % Update solution
 P = [300*u(1)^2+400*u(1)*u(2)-200*u(2)^2+150*u(1)-100*u(2);  % Nonlinear func.
      200*u(1)^2-400*u(1)*u(2)+200*u(2)^2-100*u(1)+100*u(2)];
 R = f - P;                                                          % Residual
 conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);                         % Convergence
 c = abs(0.9-u(2))/abs(0.9-uold(2))^2;                   % Convergence constant
 uold = u;                                  % Store solution for next iteration
 iter = iter + 1;                                   % Increse iteration counter
 fprintf('\n %3d  %7.5f %7.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
end                                      % End of Newton-Raphson iteration loop
%
% Example 2.6 Nonlinear algebraic equation (secant method)
tol = 1.0e-5;    iter = 0;   c = 0;       % Tolerance and initialize parameters
u = 2.0;         uold = u;                                   % Initial solution
P = u+atan(5*u); Pold = P;                                 % Nonlinear function
R = -P;          conv= R^2;                          % Residual and convergence
fprintf('\n iter       u         conv       c');                   % Print head
fprintf('\n %3d  %7.5f %12.3e %7.5f',iter,u,conv,c);
Ks = 1+5*(cos(atan(5*u)))^2;                                         % Jacobian
while conv > tol && iter < 20                       % Loop for secant iteration
 delu = R/Ks;                                    % Solve for solution increment
 u = uold + delu;                                             % Update solution
 P = u+atan(5*u);                                          % Nonlinear function
 R = -P; conv= R^2;                                  % Residual and convergence
 c = abs(u)/abs(uold)^2;                                 % Convergence constant
 Ks = (P - Pold)/(u - uold);                                           % Secant 
 uold = u; Pold = P;                        % Store solution for next iteration
 iter = iter + 1;                                   % Increse iteration counter
 fprintf('\n %3d  %7.5f %12.3e %7.5f',iter,u,conv,c);
end                                              % End of secant iteration loop
%
% Example 2.7 Two nonlinear springs (Secant method)
tol = 1.0e-5;    iter = 0;   c = 0;       % Tolerance and initialize parameters
u = [0.1; 0.1];  uold = u;                                   % Initial solution
P = [300*u(1)^2+400*u(1)*u(2)-200*u(2)^2+150*u(1)-100*u(2) % Nonlinear function
     200*u(1)^2-400*u(1)*u(2)+200*u(2)^2-100*u(1)+100*u(2)];
f = [0; 100]; R = P - f;    Rold = R;                      % Force and residual
conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);                          % Convergence
fprintf('\n iter      u1      u2         conv       c');           % Print head
fprintf('\n %3d  %7.5f %7.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
Ks = [600*u(1)+400*u(2)+150  -400*u(2)+400*u(1)-100           % Jacobian matrix
      400*u(1)-400*u(2)-100  400*u(2)-400*u(1)+100];
while conv > tol && iter < 20                       % Loop for secant iteration
 delu = -Ks\R;                                   % Solve for solution increment
 u = uold + delu;                                             % Update solution
 P = [300*u(1)^2+400*u(1)*u(2)-200*u(2)^2+150*u(1)-100*u(2);  % Nonlinear func.
      200*u(1)^2-400*u(1)*u(2)+200*u(2)^2-100*u(1)+100*u(2)];
 R = P - f;                                                          % Residual
 conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);                         % Convergence
 c = abs(0.9-u(2))/abs(0.9-uold(2))^2;                   % Convergence constant
 delR = R - Rold;                                          % Change in residual
 Ks = Ks + (delR-Ks*delu)*delu'/norm(delu)^2;            % Update secant matrix
 uold = u;  Rold = R;                       % Store solution for next iteration
 iter = iter + 1;                                   % Increse iteration counter
 fprintf('\n %3d  %7.5f %7.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
end                                              % End of secant iteration loop
%
% Example 2.8 Displacement controlled procedure
tol=1.0e-5; conv=0; u1=0; u1old=u1;       % Tolerance and initialize parameters
fprintf('\n step      u1      u2      F');                         % Print head
for i=1:9                                     % Loop for displacement increment
 u2 = 0.1*i;                                            % Increase displacement
 P = 300*u1^2+400*u1*u2-200*u2^2+150*u1-100*u2;            % Nonlinear function
 R = -P; conv = R^2;                                 % Residual and convergence
 iter = 0;                                      % Convergence iteration counter
 while conv > tol && iter < 20              % Loop for Newton-Raphwon iteration
  Kt = 600*u1+400*u2+150;                                            % Jacobian
  delu1 = R/Kt;                                  % Solve for solution increment
  u1 = u1old + delu1;                                         % Update solution
  P = 300*u1^2+400*u1*u2-200*u2^2+150*u1-100*u2;           % Nonlinear function
  R = -P; conv= R^2;                                 % Residual and convergence
  u1old = u1;                               % Store solution for next iteration
  iter = iter + 1;                                  % Increse iteration counter
 end                                     % End of Newton-Raphson iteration loop
 F = 200*u1^2-400*u1*u2+200*u2^2-100*u1+100*u2;      % Calculate reaction force
 fprintf('\n %3d  %7.5f %7.5f %7.3f',i,u1,u2,F);
end                                        % End of displacement increment loop
%
% Extension of single element example
XYZ=[0 0 0;1 0 0;1 1 0;0 1 0;0 0 1;1 0 1;1 1 1;0 1 1]; % Nodal coordinates (mm)
LE=[1 2 3 4 5 6 7 8];                                    % Element connectivity
EXTFORCE=[5 3 10; 6 3 10; 7 3 10; 8 3 10]; % External forces [Node, DOF, Value]
SDISPT=[1 1 0;1 2 0;1 3 0;2 2 0;2 3 0;3 3 0;4 1 0;4 3 0]; % Fixed displacements
TIMS=[0 0.5 0.1 0 0.5;0.5 1 0.1 0.5 1]';      % Load increments [t0 T Dt f0 fL]
MID=1; PROP=[206.9, 0.29, 0.0, 5., 35.0];      % Elastoplastic [E NU BETA H Y0]
PARAM=struct('ITRA',20,'ATOL',1.0E5,'NTOL',5,'TOL',1E-6,'TIMS',TIMS);
clc; NOUT=fopen('output.txt','w');                           % Open output file
NLFEA(PARAM,NOUT,MID,PROP,EXTFORCE,SDISPT,XYZ,LE,[],[]);        % Calling NLFEA
fclose(NOUT);                                               % Close output file
%
% Example 2.10 Two-element example
XYZ=[0 0 0; 1 0 0; 1 1 0; 0 1 0;                            % Nodal coordinates
     0 0 1; 1 0 1; 1 1 1; 0 1 1;
     0 0 2; 1 0 2; 1 1 2; 0 1 2]*0.01;
LE=[1 2 3 4 5 6 7 8; 5 6 7 8 9 10 11 12];                % Element connectivity
EXTFORCE=[9 3 1E4; 10 3 1E4; 11 3 1E4; 12 3 1E4];             % External forces
SDISPT=[1 1 0;1 2 0;1 3 0;2 2 0;2 3 0;3 3 0;4 1 0;4 3 0]; % Fixed displacements
TIMS=[0 0.8 0.4 0 0.8; 0.8 1.1 0.1 0.8 1.1]'; % Load increments [t0 T Dt f0 fL]
MID=1; PROP=[206.9E9 0.29 0 1E8 4E8]; % Elastoplastic material [E NU BETA H Y0]
PARAM=struct('ITRA',70,'ATOL',1.0E5,'NTOL',6,'TOL',1E-6,'TIMS',TIMS);
clc; NOUT=fopen('output.txt','w');                           % Open output file
NLFEA(PARAM,NOUT,MID,PROP,EXTFORCE,SDISPT,XYZ,LE,[],[]);        % Calling NLFEA
fclose(NOUT);                                               % Close output file
%
% Example 2.11 Bending of a beam (3D solid elements)
W=0.1; H=0.1; L=2; F=250;             % Dimensions of the beam and applied load
XYZ=[W 0     0;W L/5   0;0 L/5   0;0 0     0;               % Nodal coordinates
     W 0     H;W L/5   H;0 L/5   H;0 0     H;                                 %
     W 2*L/5 0;0 2*L/5 0;W 2*L/5 H;0 2*L/5 H;                                 %
     W 3*L/5 0;0 3*L/5 0;W 3*L/5 H;0 3*L/5 H;                                 %
     W 4*L/5 0;0 4*L/5 0;W 4*L/5 H;0 4*L/5 H;                                 %
     W L     0;0 L     0;W L     H;0 L     H];                                %
LE=[ 1  2  3  4  5  6  7  8;  2  9 10  3  6 11 12  7;    % Element connectivity
     9 13 14 10 11 15 16 12; 13 17 18 14 15 19 20 16;                         %
    17 21 22 18 19 23 24 20];                                                 %
FORCE=[21 3 F; 22 3 F; 23 3 F; 24 3 F];       % Applied forces [Node,DOF,Value]
SDISPT=[1 1 0;1 2 0;1 3 0;4 1 0;4 2 0;4 3 0;         % Prescribed displacements
        5 1 0;5 2 0;5 3 0;8 1 0;8 2 0;8 3 0];%               [Node, DOF, Value]
TIMS=[0.0 1.0 1.0 0 1.0]';                   % One-step linear elastic analysis
MID=0; PROP=[200E9 0.3];                              % Linear elastic material
PARAM=struct('ITRA',20,'ATOL',1.0E5,'NTOL',5,'TOL',1E-6,'TIMS',TIMS);         %
clc; NOUT = fopen('output.txt','w');                              % Output file
NLFEA(PARAM,NOUT,MID,PROP,FORCE,SDISPT,XYZ,LE,[],[]);           % Calling NLFEA
fclose(NOUT);                                             % Closing output file
[GMAX,GMIN]=PLOTELEM(XYZ,LE,'output.txt',1.0,-1,[])  % Plot undeformed geometry
[DMAX,DMIN]=PLOTELEM(XYZ,LE,'output.txt',1.0,0,[])     % Plot deformed geometry
[SMAX,SMIN]=PLOTELEM(XYZ,LE,'output.txt',1.0,2,[])    % Plot stress S22 contour
%
% P2.2 Newton-Raphson method
tol = 1.0e-5;  iter = 0;  c = 0;
u = [1; 5];    uold = u;
f = [3; 9];
P = [u(1)+u(2); u(1)^2+u(2)^2];
R = f - P;
conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);
fprintf('\n iter       u1       u2         conv       c');
fprintf('\n %3d  %8.5f %8.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
while conv > tol && iter < 20
 Kt = [1 1; 2*u(1) 2*u(2)];
 delu = Kt\R;
 u = uold + delu;
 P = [u(1)+u(2); u(1)^2+u(2)^2];
 R = f - P;
 conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);
 c = abs(3-u(2))/abs(3-uold(2))^2;
 iter = iter + 1;
 uold = u;
 fprintf('\n %3d  %8.5f %8.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
end
%
% Problem 2.3
tol = 1.0e-5;  iter = 0;  c = 0;
u = [1; 5];    uold = u;
f = [0; 9];
P = [u(1)*u(2); u(1)^2+u(2)^2];
R = f – P;
conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);
fprintf('\n iter       u1       u2         conv       c');
fprintf('\n %3d  %8.5f %8.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
while conv > tol && iter < 20
 Kt = [u(2) u(1); 2*u(1) 2*u(2)];
 delu = Kt\R;
 u = uold + delu;
 P = [u(1)*u(2); u(1)^2+u(2)^2];
 R = f – P;
 conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);
 c = abs(3-u(2))/abs(3-uold(2))^2;
 iter = iter + 1;
 uold = u;
 fprintf('\n %3d  %8.5f %8.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
end
%
% P2.4 Modified Newton-Raphson method
%
tol = 1.0e-5;  iter = 0;  c = 0;
u = [1; 5];    uold = u;
f = [3; 9];
P = [u(1)+u(2); u(1)^2+u(2)^2];
R = f - P;
conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);
Kt = [1 1; 2*u(1) 2*u(2)];
fprintf('\n iter       u1       u2         conv       c');
fprintf('\n %3d  %8.5f %8.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
while conv > tol && iter < 20
delu = Kt\R;
 u = uold + delu;
 P = [u(1)+u(2); u(1)^2+u(2)^2];
 R = f - P;
 conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);
 c = abs(3-u(2))/abs(3-uold(2))^2;
 iter = iter + 1;
 uold = u;
 fprintf('\n %3d  %8.5f %8.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
end
fprintf('\n');
%
% P2.5 Modified Newton-Raphson method
%
tol = 1.0e-5;
iter = 0;
u = [1; 5];
f = [0; 9];
P = [u(1)*u(2); u(1)^2+u(2)^2];
R = f - P;
fprintf('\n iter = %2d u = %7.3f %7.3f R= %12.3e %12.3e',...
         iter,u(1),u(2),R(1),R(2));
Kt = [u(2) u(1); 2*u(1) 2*u(2)];
while norm(R) > tol
delu = Kt\R;
u = u + delu;
P = [u(1)*u(2); u(1)^2+u(2)^2];
R = f - P;
iter = iter + 1;
fprintf('\n iter = %2d u = %7.3f %7.3f R = %12.3e %12.3e',...
         iter,u(1),u(2),R(1),R(2));
end
fprintf('\n');
%
% P2.6 Broyden method
%
tol = 1.0e-7;  iter = 0;  c = 0;
u = [1; 5];    uold = u;
f = [3; 9];
P = [u(1)+u(2); u(1)^2+u(2)^2];
R = P - f;     Rold = R;
conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);
Kt = [1 1; 2*u(1) 2*u(2)];
fprintf('\n iter       u1       u2         conv       c');
fprintf('\n %3d  %8.5f %8.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
while conv > tol && iter < 20
 delu = -Kt\R;
 u = uold + delu;
 P = [u(1)+u(2); u(1)^2+u(2)^2];
 R = P - f;
 conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);
 c = abs(3-u(2))/abs(3-uold(2))^2;
 delR = R - Rold;
 Kt = Kt + (delR-Kt*delu)*delu'/norm(delu)^2;
 uold = u;  Rold = R;
 iter = iter + 1;
 fprintf('\n %3d  %8.5f %8.5f %12.3e %7.5f',iter,u(1),u(2),conv,c);
end
fprintf('\n');
%
% P2.7 Incremental force method
%
tol = 1.0e-5;   u = [1; 5];
fprintf('\n  inc       u1       u2       f1       f2');
for i=1:5
 f = i*0.2*[3; 9];
 P = [u(1)+u(2); u(1)^2+u(2)^2];
 R = f - P;
 conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);
 iter = 0;
 while norm(R) > tol && iter < 20
  Kt = [1 1; 2*u(1) 2*u(2)];
  delu = Kt\R;
  u = u + delu;
  P = [u(1)+u(2); u(1)^2+u(2)^2];
  R = f - P;
  conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);
  iter = iter + 1;
 end
 fprintf('\n %3d  %8.5f %8.5f %8.5f %8.5f',i,u(1),u(2),f(1),f(2));
end
fprintf('\n');
%
% P2.8 Incremental force method
%
tol = 1.0e-5;
iter = 0;
u = [1; 5];
f0 = [0; 9];
for I = 1:5
 f = I*0.2*f0;
 P = [u(1)*u(2); u(1)^2+u(2)^2];
 R = f - P;
 fprintf('\n iter = %2d I = %2d u = %7.3f %7.3f R = %12.3e %12.3e',...
          iter,I,u(1),u(2),
 R(1), R(2));
 while norm(R) > tol
  Kt = [u(2) u(1); 2*u(1) 2*u(2)];
  delu = Kt\R;
  u = u + delu;
  P = [u(1)*u(2); u(1)^2+u(2)^2];
  R = f - P;
  iter = iter + 1;
  fprintf('\n iter = %2d I = %2d u = %7.3f %7.3f R = %12.3e %12.3e',...
           iter,I,u(1),u(2),
  R(1),R(2));
 end
end
fprintf('\n');
%
% P2.9 Nonlinear strain bar
%
fprintf('\n inc        F       u  strain');
tol = 1.e-5;  u = 0;
for i=1:10
 strain = 0.5*u^2 + u;
 f = 0.1*i;
 iter = 0;
 P = u^3+3*u^2+2*u;
 R = f - P;
 conv= R^2/(1+f^2);
 while conv > tol && iter < 20
  Kt = 3*u^2+6*u+2;
  delu = R/Kt;
  u = u + delu;
  strain = 0.5*u^2 + u;
  P = u^3+3*u^2+2*u;
  R = f - P;
  conv= R^2/(1+f^2);
  iter = iter + 1;
 end
 fprintf('\n %3d  %7.5f %7.5f %7.5f',i,f,u,strain);
end
%
% P2.10 Nonlinear strain bar (secant method)
%
tol = 1.0e-5;      iter = 0;
u = 0.0;           uold = u;   f = 1;
P = u^3+3*u^2+2*u; Pold = P;
R = f - P;         conv= R^2/(1+f^2);
strain = 0.5*u^2 + u;
fprintf('\n iter       u  strain         conv');
fprintf('\n %3d  %7.5f %7.5f %12.3e',iter,u,strain,conv);
Ks = 3*u^2+6*u+2;
while conv > tol && iter < 20
 delu = R/Ks;
 u = uold + delu;
 strain = 0.5*u^2 + u;
 P = u^3+3*u^2+2*u;
 R = f - P;
 conv= R^2/(1+f^2);
 Ks = (P - Pold)/(u - uold);
 uold = u;
 Pold = P;
 iter = iter + 1;
 fprintf('\n %3d  %7.5f %7.5f %12.3e',iter,u,strain,conv);
end
%
% P2.11 Elasto-plastic bar
%
fprintf('\n inc         F       u     stress');
tol = 1.e-5;  u = 0;  uY = 0.002;
for i=1:10
 f = i*5000;
 iter = 0;
 if u < uY P = 2E7*u;
 else      P = 2E7*uY + 2E6*(u-uY); end
 R = f - P;
 conv= R^2/(1+f^2);
 while conv > tol && iter < 20
  if u < uY Kt = 2E7;
  else      Kt = 2E6; end
  delu = R/Kt;
  u = u + delu;
  if u < uY P = 2E7*u;
  else      P = 2E7*uY + 2E6*(u-uY); end
  R = f - P;
  conv= R^2/(1+f^2);
  iter = iter + 1;
 end
 if u < uY stress = 2E11*u;
 else      stress = 4E8 + 2E10*(u-uY); end
 fprintf('\n %3d  %8.1f %7.5f %10.3e',i,f,u,stress);
end
%
% P2.12 Three nonlinear springs
%
tol = 1.0e-5;  iter = 0;
u = [0; 0];    uold = u;
f = [0; 100];
P = [50*u(1)^2+200*u(1)*u(2)-100*u(2)^2+1200*u(1)-500*u(2)
     100*u(1)^2-200*u(1)*u(2)+100*u(2)^2-500*u(1)+500*u(2)];
R = f - P;
conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);
fprintf('\n iter      u1      u2         conv');
fprintf('\n %3d  %7.5f %7.5f %12.3e',iter,u(1),u(2),conv);
while conv > tol && iter < 20
 Kt = [100*u(1)+200*u(2)+1200   200*u(1)-200*u(2)-500
       200*u(1)-200*u(2)-500   -200*u(1)+200*u(1)+500];
 delu = Kt\R;
 u = uold + delu;
 P = [50*u(1)^2+200*u(1)*u(2)-100*u(2)^2+1200*u(1)-500*u(2)
      100*u(1)^2-200*u(1)*u(2)+100*u(2)^2-500*u(1)+500*u(2)];
 R = f - P;
 conv= (R(1)^2+R(2)^2)/(1+f(1)^2+f(2)^2);
 uold = u;
 iter = iter + 1;
 fprintf('\n %3d  %7.5f %7.5f %12.3e',iter,u(1),u(2),conv);
end
%
% P2.13 Nonlinear modulus bar
%
fprintf('\n inc iter       F       u       stress');
tol = 1.e-5;  u = 0;
for i=1:10
 f = i*2500;
 iter = 0;
 stress = 1E9*(1-u)*u;
 P = stress*1E-4;
 R = f - P;
 conv= R^2/(1+f^2);
 while conv > tol && iter < 20
  Kt = 1E5*(1-2*u);
  delu = R/Kt;
  u = u + delu;
  stress = 1E9*(1-u)*u;
  P = stress*1E-4;
  R = f - P;
  conv= R^2/(1+f^2);
  iter = iter + 1;
 end
 fprintf('\n %3d %3d  %7.1f %7.5f %12.3e',i,iter,f,u,stress);
end
%
% Problem 2.15 Simple shear of a square block
XYZ=[0 0;1 0;0 1;1 1];                                      % Nodal coordinates
LE=[1 2 4 3];                                            % Element connectivity
EXTFORCE=[3 1 5000;4 1 5000];                                 % External forces
SDISPT=[1 1 0;1 2 0;2 1 0;2 2 0;3 2 0;4 2 0];                    % Fixed displ.
TIMS=[0 1 1 0 1]';                            % Load increments [t0 T Dt f0 fL]
MID=0; PROP=[200E9 0.3];                       % Linear elastic material [E NU]
PARAM=struct('ITRA',10,'ATOL',1.0E5,'NTOL',6,'TOL',1E-6,'TIMS',TIMS);         %
clc; NOUT=fopen('output.txt','w');                           % Open output file
NLFEA(PARAM,NOUT,MID,PROP,EXTFORCE,SDISPT,XYZ,LE,[],[]);        % Calling NLFEA
fclose(NOUT);                                               % Close output file
%
% Problem P2.16 2D Bending of a beam (5-elements)
XYZ=[0 0;1 0;2 0;3 0;4 0;5 0;0 1;1 1;2 1;3 1;4 1;5 1];      % Nodal coordinates
LE=[1 2 8 7;2 3 9 8;3 4 10 9;4 5 11 10;5 6 12 11];       % Element connectivity
FORCE=[6 1 10000;12 1 -10000];                                 % Applied forces
SDISPT=[1 1 0;1 2 0;7 1 0];                           % Prescribed displacement
MID=0; PROP=[2E6, 0.3];                  % Linear elastic material PROP=[E, NU]
TIMS=[0 1 1 0 1]';                                             % Load increment
PARAM=struct('ITRA',10,'ATOL',1E5,'NTOL',5,'TOL',1E-6,'TIMS',TIMS);% Parameters
clc; NOUT=fopen('output.txt','w');                                % Output file
NLFEA(PARAM,NOUT,MID,PROP,FORCE,SDISPT,XYZ,LE,[],[]);              % Call NLFEA
fclose(NOUT);                                               % Close output file
[SMAX,SMIN]=PLOTELEM(XYZ,LE,'output.txt',1.0,1,1)      % Plot deformed geometry