function [y,t]=rungekuttaImplicit(f,ta,ya,te,N,M) 
 h=(te-ta)/N;
 h6=h/6;
 t = ta:h:te; % Gitter anlegen
 y=zeros(length(ya),N+1);
 y(:,1) = ya; % Startwert schreiben

 for i=1:N % durchlaufe Gitterpunkte
    %setze Startwerte fuer Iterationsverfahren
    k1=zeros(length(ya),1);
    k2=zeros(length(ya),1);
    %starte Iterationsverfahren
    for j=1:M
      k1alt=k1;  % alte Werte merken fuer Iteration
      k2alt=k2;
      k1=f(t(i),y(:,i)+h6*(k1alt-k2alt));
      k2=f(t(i)+h/2,y(:,i)+h6*(k1alt+2*k2alt));
    end
    k3=f(t(i)+h,y(:,i)+h6*(k1+5*k2));
    y(:,i+1)=y(:,i)+h6*(k1+4*k2+k3);
 end
end