% Antti Hannukainen 2.3.2006
%
% Compute all values needed in the continuation 
% of the periodic solution
%
% i.e. the Jacobian rising from the Midpoint - discretized of the equation
%
% x'           = tau*f(x,a)
% g(x(0),x(1)) = 0
% sigma(x)     = 0
%
%
% [fu, Dfu] = per_Df(u,u0,t,f,g,t)
%
% u  = the current point = [ ux ; tau ; a]
% u0 = previous iterate 
% t  = mesh in time  = [ 0  h1  .... 1];
% 
% f = f(ux,a); 
% sigma = sigma(u,u0,t)
%


function [fu,Dfu] = per_Df(u,u0,t,f,sigma)


% convert u - into more suitable form 
N = length(u);

% Current point
U  = u( 1:(N-2));
tau  = u(N-1);
a    = u(N);

% Old point (Needed in the  phase condition)
U0 = u0(1:(N-2));
tau0 = u0(N-1);
a0   = u0(N);

% Define the boundary value function 
g = @(x0,xm)(x0-xm);

% Compute as much as possible with the existing bvp-code
%

% Dux
M = bvp_Dfu( U, @(t,x)(tau*f(x,a)), g, t);  

% Dtau
[Qtau, xxtau, fftau, bb ] = bvp_fu(U, @(t,x)(f(x,a)), g, t);

% Da 
[Qa, xxa, ffa] = bvp_fu( U, @(t,u)(Dapprox( @(x)(tau*f(u,x)),a )),g,t);

% Ds
Ds = Dapprox( @(z)(phase_condition(z,U0,t)), U);


% Combine into Jacobian and into value
Dfu = full( [M  -fftau -ffa ; Ds [0 0]] );

fu = full( [xxtau-tau*fftau+bb; phase_condition(U,U0,t)]);  

