function J=JacG(theta)
%
global T alpha beta
m=size(theta,1);
h=T/(m+1);
e=ones(m,1);
cc=cos(theta);
J=-1/(h^2)*spdiags([e -2*e+h^2*cc e],-1:1,m,m);