f = inline('x/y');
%f = inline('-y/(x+1)');
y = inline('sqrt(1+x.^2)');
%y = inline('1./(x+1)')
dy = inline('x./sqrt(1+x.^2)');
d2y = inline('1./((1+x.^2).^1.5)');
d3y = inline('-3*x./((1+x.^2).^2.5)');


h = 0.1;
T = 2;
N = T/h;

% Eulero esplicito
yeul(1) = 1;
for i=1:N
 yeul (i+1) = yeul(i) + h*f((i-1)*h,yeul(i));
end

% Crank Nicolson   
ycn(1) = 1;
for i=1:N
 ycn_0 = ycn(i);
 scarto = 1;
 iter = 0;
 while (scarto > h^3/10 && iter < 10)
   ycn_k = ycn(i) + h*(f((i-1)*h,ycn(i))+f(i*h,ycn_0))/2;
   scarto = abs(ycn_k -ycn_0);
   ycn_0 = ycn_k;
    iter = iter +1;
end
iter
ycn(i+1) = ycn_k;
end
x = (1:N)*h;
figure(1);
plot(x,yeul(2:N+1),'bo',x,ycn(2:N+1),'r*',x,y((1:N)*h),'b-')
figure(2)
semilogy(x,abs(yeul(2:N+1)-y(x)),'bo',x,abs(ycn(2:N+1)-y(x)),'r*')
errEU = h/2*d2y(0)
errCN = h^2/12*d3y(0.5)
