function y = fun1( x )
global l0 b k m g

  y = k * (sqrt( b^2 + x^2 ) - l0) * x / sqrt( b^2 + x^2 ) - m * g;