Tend = 30; u0 = [0;0];

function curr = Diode(u)
  Rd = 10; us = 0.7;
  if (u>=-us)  curr = 0;
  else curr = Rd*(u+us);
  endif
endfunction

function y = circuit(t,u)
  C1 = 1; C2 = 1;
  y  = [-1/C1*(Diode(u(1)-10*sin(t))-Diode(u(2)-u(1)));
        10*cos(t)-1/C2*Diode(u(2)-u(1))];
endfunction


t0 = cputime();
[t,u] = rk45('circuit',0,Tend,u0,1e-5,1e-5); % Runge Kutta adaptiv
timer = cputime()-t0
figure(1);
plot(t,u(:,2)'.')
grid on; xlabel('time'); ylabel('tension');

tFix = linspace(0,Tend,100);
t0 = cputime();
[tFix,uFix] = ode_Runge('circuit',tFix,u0,1); % Runge Kutta
timer = cputime()-t0
figure(2);
plot(tFix,uFix(:,2),t,u(:,2))
grid on; xlabel('time'); ylabel('tension');
legend('u fix','u adapt')