% [y,x] = cranknicholson(f,a,ya,n);
% Implicitna Crank-Nicholsonova metoda za resevanje zacetnega problema y'=f(x,y), y(x0)=y0.
% Vhod:
% f : funkcija f(x,y) -- desna stran diferencialne enacbe y'=f(x,y)
% a : interval [a, b] na katerem iscemo numericno resitev
% ya : vektor zacetnih pogojev v levem krajiscu
% n : stevilo podintervalov danega intervala
% Izhod:
% x : delitev intervala [a, b]
% y : matrika vrednosti funkcije y(x) v notranjih tockah intervala in na robu
function [y,x] = cranknicholson(f,a,ya,n,e);
  m = 20;
  if nargin < 5
    e = 1e-10;
  end;
  h = (a(2)-a(1))/n;
  y = zeros(length(ya),n+1);
  y(:,1) = ya(:);
  x = linspace(a(1),a(2),n+1);
  for i = 2:n+1
    ff = f(x(i-1),y(:,i-1));
    y0 = y(:,i-1);
    for j = 1:m
      y1 = y(:,i-1)+h/2*(ff+f(x(i),y0));
      if abs((y0-y1)/y1) < e break; end;
      y0 = y1;
    end;
    y(:,i) = y1;
  end;
