function xk1=radau3_step(f,xk,tk,h,dxguess)
  function ff=f_xx(k12,tk)
    n=length(k12)/2;
    k1a=k12(1:n);
    k2a=k12(n+1:2*n);
    ff=[f(xk+5*h/12*k1a-h/12*k2a,k1a,tk+h/3);f(xk+3*h/4*k1a+h/4*k2a,k2a,tk+h)]
  endfunction
  kprev=[dxguess;dxguess];
  J=jacobiant(f_xx,kprev,tk+h);//Jacobiano del Sistema
  invH=pinv(J);//inversa del Hessiano
  kk=kprev-invH*f_xx(kprev,tk+h);//iteracion de Newton
  l=0;
  while norm(kk-kprev)>1e-6*(norm(kprev)+1e-3)&l<20
    l=l+1;
    kprev=kk;
    J=jacobiant(f_xx,kprev,tk+h);//Jacobiano del Sistema
   invH=pinv(J);//inversa del Hessiano
   kk=kprev-invH*f_xx(kprev,tk+h);//iteracion de Newton
end  
  n=length(kk)/2;
  k1=kk(1:n);
  k2=kk(n+1:2*n);
  xk1=xk+h/4*(3*k1+k2);
endfunction

function [t,x]=radau3(f,x0,ti,tf,h)
//metodo de Backward Euler con iteracion de Newton
//exec('jacobiant.sci')
  t=[ti:h:tf];
  x=zeros(length(x0),length(t));
  x(:,1)=x0;
  xn=beulerdae_step(f,x0,t(1),h);
  for k=1:length(t)-1
   dxguess=(xn-x(:,k))/h;
   x(:,k+1)=radau3_step(f,x(:,k),t(k),h,dxguess);
   xn=x(:,k+1);
  end
endfunction


