%
%  lect03.m -- ODE script
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

global A B

A = 1.1;  B = 1.0;  tol = 1e-2;  Up0 = 0.3704;  Xode = 'Xsteady';  
%A = 4.1;  B = 1.0;  tol = 1e-2;  Up0 = 0.7320;  Xode = 'Xsteady2';  

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

%  IVs, final time, method & options

IVs = [0 Up0];  T = 10;  method = 'ode45';
options = odeset('reltol',tol*1e-3,'abstol',tol*1e-6,'refine',1, ...
		'maxstep',T/20,'stats','off','events','on');

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

%  ODE solve 

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

%  uniform grid interpolation

N = 32;  L = t(end);  dx = L/N;

xg = dx:dx:L-dx;  
ug = interp1(t,y(:,1),xg,'cubic');

xp = [0 xg L];  up = [0 ug 0];

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

clf;subplot(2,1,1)

plot(t,y(:,1),'kx');  hold on;  plot(t,y(:,1),'b')
plot(xp,up,'ko');  plot(xp,up,'r')

title('\bf steady-state (u_{ode}=blue, u_{unif}=red)')
xlabel('\bf x-axis');  ylabel('\bf u-axis');

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

%  eigenvalue analysis
%    - build matrix (help spdiags)

d2 = diag(ones(N-2,1),1) - 2*diag(ones(N-1,1),0) + diag(ones(N-2,1),-1);
Ma = d2/(dx^2) - 3*B*diag(ug.^2,0);

%  eigenfunction plot

[Vec Evs] = eig(Ma);
[stab index] = max(diag(Evs));
Efun = Vec(:,index)/max(Vec(:,index));
disp(['stability eigenvalue = ' num2str(stab+A)]);  disp('  ')

subplot(2,1,2)
plot(xp,[0 Efun' 0],'g');  hold on;  plot(xp,[0 Efun' 0],'ko')

title(['\bf stability eigenfunction'])
xlabel('\bf x-axis');  ylabel('\bf w(x)');
axis([0 Tend 0 1.1])
text(1,0.2,['eigenvalue = ' num2str(stab+A)])