% Antti Hannukainen 1.3
%
% In BVP - solution process one needs to compute Jacobian DFu of
% the discretized equation. This function computes and
% Jacobian and it is ment to be used as
%
% Dfu = @(u)( bvp_Dfu(u,f,g,t))
%
% Where u is the point, where the Jacobian is to be computed, f and
% g are function handles to BVP - funcions f and g. The needed
% jacobians are computed with central difference approxiation
% (Dapprox)
%
%  t = ia a vector of discrete time steps
%
%
%  u = [ x0  x1  x2  x3  ... xm];
%
%

function DFu = bvp_Dfu(u,f,g,t)

x = bvp_u2x(u,length(t));

% Define helper - functions to simplify computation of the Jacobian. 
Df  = @(t,x)(Dapprox( @(z)(f(t,z)), x));
Dgx = @(x0,xm)( Dapprox( @(z)( g(z , xm) ), x0) );
Dgy = @(x0,xm)( Dapprox( @(z)( g(x0 , z) ), xm) );

% Define some constants
Dim = length( f(t(1),x(:,1)) );
N = length(t);

M = sparse(Dim*N);

for i=1:1:(N-1)
  
  dt = t(i+1) - t(i);
  
  abu1 = [(Dim*(i-1)+1):( Dim*i)];
  abu2 = [(Dim*i+1):(Dim*(i+1) )];
  
  % ADD HERE CODE TO COMPUTE DFU
  
end

% ADD CODE TO Set Boundary Conditions


abu1 = [(Dim*(N-1)+1):(Dim*N)];
abu2 = [1:Dim];
  
M(abu1,abu2) = B0;

abu2 = [(Dim*(N-1)+1):(Dim*N)];

M(abu1,abu2) = B1;

DFu = M;
