clear;
format long;

ns=400;
n=linspace(1,ns,ns);
cc=0; h=sqrt(1/cc);
%u1=4*atan(exp(h*(n-ns/2)));
%u1=4*atan(exp(h*(n-ns/2-0.5)));

u1=zeros(1,ns);
for i=ns/2+1:ns
  u1(i)=2*pi;
end;
% include line below for the unstable solution.
u1(ns/2)=pi;


it=1;

while (cc<4)

u=zeros(1,ns); ee=ones(ns,1);

while (norm(u-u1)>1e-08)

u=u1;

f=cc*[0,diff(u,2),0]-sin(u);

jj=spdiags([cc*ee,-2*cc*ee-cos(u'),cc*ee],-1:1,ns,ns);

du=(-jj\f')';

norm(du);

u1=u+du;

end;

jj1=full(jj);
[vv,dd]=eig(jj1);
dd1=diag(dd);

c_all(it)=cc;
u_all(:,it)=u;
d_all(:,it)=dd1;
it=it+1;
cc=cc+0.01
end;

fid=fopen('sg_c05.dat','w')
for i=1:ns
fprintf(fid,'%6.10f\n',u_all(i,51));
end;