function bessel_cn

h=0.01;
n=0;

t=0:h:10;
[l,m]=size(t);

x(:,1)=[0;1];
A=[1, -h*(n/t(2)-t(2));-h/t(2),1];
x(:,2)=inv(A)*x(:,1);


for k=3:m
    A=[1, -h*(n/t(k)-t(k))/2;-h/(2*t(k)),1];
    B=[1, h*(n/t(k)-t(k))/2;h/(2*t(k)),1];
    x(:,k)=inv(A)*B*x(:,k-1);
end

figure(2)
plot(t,x(2,:),'b',t,bessel(n,t),'r')