%
%  lect02.m -- ODE script
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

global A B
A = 1.0;  B = 1.0;  tol = 1e-2;

disp(' ')
Up0 = input('initial derivative (Up0):  ');

if (abs(Up0) > sqrt(0.5*A^2/B)) 
    disp('first integral error');  disp(' ')
    break
end

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

Xode = 'Xsteady';  IVs = [0 Up0];  T = 10;  method = 'ode45';

options = odeset('reltol',tol*1e-3,'abstol',tol*1e-6,'refine',1, ...
		'maxstep',T/20,'stats','on','events','on');

disp(' ');  disp([method ' for ' Xode]);  disp(' ')

%  ODE solve 

[t y] = feval(method,Xode,[0 T],IVs,options);

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

clf;subplot(2,1,1)

plot(t,y(:,1),'kx')
hold on
plot(t,y(:,1),'b')

plot(t,y(:,2),'kx')
plot(t,y(:,2),'r')

title('\bf steady-state (u=blue, v=red)')
xlabel('\bf x-axis');  ylabel('\bf u,v-axis');

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

subplot(2,1,2)
plot(t(2:end),diff(t),'kx')
title(['\bf timesteps for ' method])
xlabel('\bf x-axis');  ylabel('\bf \Delta x');
axis([0 Tend 0 1.1*max(diff(t))])

%  first integral

check = (y(:,2).^2)/2 + A*(y(:,1).^2)/2 - B*(y(:,1).^4)/4 ;
disp(' ')
disp(['1st-integral variation = ' num2str(max(check)-min(check))])

data = [data; Up0 Tend]