function DF = brachijac ( x );
% BrachiJac
%
% Jacobian for trapezoidal rule applied to Brachistochrone problem
%
% Dieses m-file ist Teil der Arbeitsmaterialien des 
%   Schuelerseminars "Optimale Bahnplanung von Robotern"
% am Fachbereich Mathematik und Informatik der
% Martin-Luther-Universitaet Halle-Wittenberg, Deutschland
%
% Author :      Prof. Dr. M. Arnold, martin.arnold@mathematik.uni-halle.de
% Version of :  Feb 8, 2006

% -> boundary data
y_0 =  0;
y_B = -0.8;

u_0 = (-y_0)^(3/2);
u_B = (-y_B)^(3/2);

% -> read input
N = length(x);
h = 1/N;

u = [ u_0; x(1:(N-1)); u_B ];
C = x(N);

% -> allocate memory
DF = zeros ( N, N );

% -> evaluate elements of Jacobian
u13 = zeros ( size(u) );
u13(u>=0) = u(u>=0).^(1/3);
u13(u<0)  = -(-u(u<0)).^(1/3);

DF(1,1) = 1/h + 1/(4*sqrt(C^2-(u(1+1)^2)^(1/3))*u13(1+1));
for i=2:N-1,
    DF(i,i-1) = -1/h + 1/(4*sqrt(C^2-(u(1+i-1)^2)^(1/3))*u13(1+i-1));
    DF(i,i)   =  1/h + 1/(4*sqrt(C^2-(u(1+i)^2)^(1/3))  *u13(1+i));
end;
DF(N,N-1) = -1/h + 1/(4*sqrt(C^2-(u(1+N-1)^2)^(1/3))*u13(1+N-1));

for i=1:N,
    DF(i,N) = - ( 3/2*C/sqrt(C^2-(u(1+i-1)^2)^(1/3)) + ...
                  3/2*C/sqrt(C^2-(u(1+i  )^2)^(1/3)) ) / 2;
end;

