% Antti Hannukainen 
%
% 

function [Fu,xx,ff,bb] = bvp_fu(u,f,g,t)

% Convert x to simpler format
x = bvp_u2x(u,length(t));

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


for i=1:1:(N-1)
  
  dt = t(i+1) - t(i);
  
  abu1 = [(Dim*(i-1)+1):( Dim*i)];
  
  ff(abu1,1)   =  dt*f( t(i) + dt/2 ,0.5*(x(:,i)+x(:,i+1)) );
  xx(abu1,1)   =  x(:,i+1) - x(:,i);

  % ADD FU HERE !!!
  
end

% Set Boundary Conditions
abu1 = [(Dim*(N-1)+1):(Dim*N)];
b  =  g( x(:,1), x(:,N));

Q(abu1,1)=b;


bb(abu1,1) = b;
ff(abu1,1) = 0;
xx(abu1,1) = 0;

Fu=Q;

