global l0 b k m g

echo off

l0 = 5.0;
b = 5.0;
m = 1.0;
g = 10.0;

% Numerical solution.
clf
view(3)
hold on;
tic
for k = 1:0.1:10
   for l0 = 1:2:20
      xz = fzero( 'fun1', 0 );
      pull = (sqrt( b^2 + xz^2 ) > l0);
      color = pull * 'gv' + (~pull) *  'c^';
      plot3( k, l0, xz, char( color ) )
  end
end
toc
hx = xlabel( 'k [N/m]' );
set( hx, 'fontsize', 14 );
set( hx, 'color', 'r' );

hy = ylabel( 'l0 [m]' );
set( hy, 'fontsize', 14 );
set( hy, 'color', 'r' );

hz = zlabel( 'x [m]' );
set( hz, 'fontsize', 14 );
set( hz, 'color', 'r' );


h = title( 'Equilibrium position - k, l0 dependence');

echo on
get( h )
set( h, 'fontsize', 14 );
set( h, 'color', 'b' );

echo off
% $$$ % Symbolic solution.
% $$$ fstring = 'k * (sqrt( b^2 + x^2 ) - l0) * x / sqrt( b^2 + x^2 ) - m * g';
% $$$ 
% $$$ fsym = sym( fstring );
% $$$ 
% $$$ x = solve( fsym, 'x' );
% $$$ x = simple( x )
% $$$ xx = subs( x )
