scf(0);
clf()
t=0:0.001:20;

y = besselj(0,t);
plot(t,y,'k-')

n = size(y,2);
y2=y(2:n);y2(n)=0;
y3 = y .* y2;
z = find(y3<0);

scf(1);
clf()
s=0:0.01:1;
c=['k:','k-','b:','b-','r:','r-','g:','g-'];
for j = 1:size(z,2);
	plot(s,besselj(0,t(z(j))*s),c(j));
end;

jnorm = zeros(size(z,2));
for j = 1:size(z,2);
	jnorm(j) = inttrap(s,besselj(0,t(z(j))*s).^2 .*s );
end;


jcoef = zeros(size(z,2));
for j = 1:size(z,2)
	jcoef(j) = inttrap(s,besselj(0,t(z(j))*s).*(1-s.^2).*s);
end;

jc=jcoef./(jnorm)

scf(2);
clf()
plot(s,1-s.^2,'rx-');
plot(s,besselj(0,t(z(1))*s)*jc(1),'c:');
plot(s,besselj(0,t(z(1))*s)*jc(1)+besselj(0,t(z(2))*s)*jc(2),'c-')
plot(s,besselj(0,t(z(1))*s)*jc(1)+besselj(0,t(z(2))*s)*jc(2)+besselj(0,t(z(1))*s)*jc(3),'g:')
plot(s,besselj(0,t(z(1))*s)*jc(1)+besselj(0,t(z(2))*s)*jc(2)+besselj(0,t(z(1))*s)*jc(3)+besselj(0,t(z(1))*s)*jc(4),'ko-')

t(z)