%
%  code29.m -- gumdrop oscillator (26 nov 2001 - djm)
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

disp(' ')
Up = input('initial derivative (U''):  ');

%  ODE file, IVs, final time, method & options

Xode = 'code29a';  IVs = [0 Up];  T = 100;  method = 'ode45';  
tol = 1e-12;

options = odeset('reltol',tol,'abstol',tol,'refine',1, ...
		'maxstep',T/20,'stats','off','events','on');

%  forward & backward ODE solve

[tf yf] = feval(method,Xode,[0  T],IVs,options);
[tb yb] = feval(method,Xode,[0 -T],IVs,options);

y = [yb(end:-1:2,:) ; yf];  t = [tb(end:-1:2) ; tf];
Tend = tf(end);  period = tf(end)-tb(end);

%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%
%  solution plot

figure(2);  clf;
subplot(2,1,1)

plot(t,y(:,1),'b');  hold on
plot(t,y(:,2),'r')

title('\bf gumdrop oscillation (ode45)')
xlabel('\bf t-axis');  ylabel('\bf u(t) in blue, u\prime(t) in red');

amax = max(max(abs(y)));
axis([-Tend Tend -1.1*amax 1.1*amax])

%  first integral

check  = (y(:,1).^2) + y(:,2) - log(2*y(:,2)+1)/2 ;
check0 = (IVs(1).^2) + IVs(2) - log(2*IVs(2)+1)/2 ;
disp(' ')
err = max(abs(check-check0));
disp(['1st-integral variation = ' num2str(err,8)]);
disp(['period = ' num2str(period,8)]);
disp(' ');

subplot(2,1,2)

if(0==1)
  plot(t,check-check0)
  axis([-Tend Tend -1.1*err 1.1*err])
  title('\bf first integral error')
  xlabel('\bf t-axis');  ylabel('\bf error');
else
  plot(t,y(:,1),'b');  hold on	
  plot(t,2*y(:,2).*y(:,1),'r')
  plot(t,-(2*y(:,2)+1).*y(:,1),'g');
  amax = max(abs((2*y(:,2)+1)).*y(:,1));
  axis([-Tend Tend -1.1*amax 1.1*amax])
  title('\bf balance analysis (ode45)')
  xlabel('\bf t-axis');  ylabel('\bf u=b, 2uu\prime =r, u\prime\prime = g');
end

figure(1);
subplot(1,2,2)
x = -5:0.1:5;  plot(x,Up - x.^2,'g');  hold on
plot(y(:,1),y(:,2),'r')

title('\bf gumdrop phase plane (ode45)')
xlabel('\bf u(t)');  ylabel('\bf u\prime(t)');
axis([min(y(:,1)) max(y(:,1)) -0.5 1.1*Up])
axis equal
hold off