Displacements not correct issue

[hlyr]

Hello,

I am testing the simple problem of a cantilever with a concentrated load at its tip, using various different beam elements. For some reason the deflection at the tip starts to deviate from the linear elastic one for the most cases of varying the load value. At the end I print the data info for the node at the tip, at the same time, print the the value of the tip deflection due to linear elastic response ( P*L^3/(3EI) ) - along with the shear component, although this is negligible. The value from the built-in print command and the one I typed to print should coincide, but thats only the case for certain values of loading P (1.0 kN for example).

Cantilever is horizontal, of length 0.5m(eters) with rectangular cross-section of dimension d=b=2cm. By varying Py I get different deflections from print, compared to the expected linear elastic response.
Below is the code I used:

#########################################################
wipe
model BasicBuilder -ndm 2 -ndf 3

# Use elastic beam column element
set ele elasticBeamColumn

# Units
set m 1.
set GPa 10e6
set cm [expr $m/100]
set kN 1.

# Mechanical and geometric parameters
set L [expr 0.5*$m]
set E [expr 200*$GPa]
set ni 0.3; #Poisson's ratio
set G [expr $E/(2.0*(1 + $ni))]
set d [expr 2.0*$cm]
set b $d
set A [expr $d*$b]
set I [expr $b*$d**3/12.0]

set steps 100; #number of steps

# Tip Loads
set Px [expr 0.0*$kN]
set Py [expr 5.0*$kN]
set dP [expr $Py/$steps]

# nodes
node 1 0.0 0.0
node 2 0.1 0.0
node 3 0.2 0.0
node 4 0.3 0.0
node 5 0.4 0.0
node 6 0.5 0.0

fix 1 1 1 1

geomTransf Linear 1

# Five elements
element $ele 1 1 2 $A $E $I 1
element $ele 2 2 3 $A $E $I 1
element $ele 3 3 4 $A $E $I 1
element $ele 4 4 5 $A $E $I 1
element $ele 5 5 6 $A $E $I 1

timeSeries Linear 1

pattern Plain 1 1 {
load 6 $Px $Py 0.0
}

# Solve
system BandGeneral; # system of equations
numberer Plain; # enumerate DOFs
constraints Plain
test NormDispIncr 1.0e-6 20
integrator LoadControl $dP
algorithm Newton
analysis Static

analyze $steps

print node 6

# Print linear elastic deflection at tip
set f_disp [expr $Py*$L**3/(3*$E*$I)]
set s_disp [expr (5.0/6.0)*$Py*$L/($G*$A)]

puts "Displacement at tip is:"
puts "Due to flexure: $f_disp m"
puts "Due to shear: $s_disp m"
puts "Total displacement is: [expr $f_disp + $s_disp] m"