% 23.04.2001, c
fig = 1;
figure( fig );
clf

x0 = [0, 0]';
v = 1;

xa0 = [10, 0]';

options = odeset( 'OutputFcn', 'odeplot' );

%fnXa = '0';
fnXa = '-t/3';
%fnYa = '5 * sin(t)';
fnYa = '4 * sin(t) + 3 * cos(2*t)';
tspan = [0 20];
%[tt, xx] = ode23( 'psi_rhs', tspan, x0, options, xa0, v, fnXa, fnYa );
[tt, xx] = ode45( 'psi_rhs', tspan, x0, options, xa0, v, fnXa, fnYa );
%[tt, xx] = ode113( 'psi_rhs', tspan, x0, options, xa0, v, fnXa, fnYa );

pause

clf

l = norm( x0 - xa0 );
axis( [x0(1), xa0(1), x0(2) - l / 2, x0(2) + l / 2] );
hold on;
hx = plot( x0(1), x0(2), 'b*', 'erasemode', 'xor' );
hxa =  plot( xa0(1), xa0(2), 'r+', 'erasemode', 'xor' );
options = [];

dt = 0.05;
jj = 1;
jj0 = 1;
for t = tspan(1):dt:tspan(2);
  while (tt(jj) < t)
    jj0 = jj;
    jj = jj + 1;
  end
  if (jj == jj0)
    x = xx(jj,:);
  else
  x = ((t - tt(jj0)) * xx(jj,:) + (tt(jj) - t) * xx(jj0,:))...
      / (tt(jj) - tt(jj0));
  end
% $$$   [jj0 jj -1 t tt(jj) tt(jj0); xx(jj,:) xx(jj0,:) x]
% $$$   [ (t - tt(jj0)), (tt(jj) - t)]
% $$$   pause
  xa = xa0 + [eval( fnXa ); eval( fnYa ) ];
  set( hx, 'marker', '.' );
  hx = plot( x(1), x(2), 'b*', 'erasemode', 'xor' );
  set( hxa, 'xdata', xa(1), 'ydata', xa(2) );
  drawnow
  tic;
  while ((toc) < 0.01)
    ;
  end
end
plot( xx(:,1), xx(:,2) );
