I've used ODE45 to solve an equation, now: How to use Runge-Kutta 4 and interpolation to find two specific values?

Hi, first of all pardon my writing (English is not my native language) and the very long post I am about to put here. But I'm seriously stuck and need help :) I have solved the following equation using ODE45, with the purpose to calculate and plot power curves (see further below):

du = @(t,u)[ u(2); (2*u(1)/(1+u(1)^2))*u(2)^2-(1+u(1)^2)*u(1)/(L0*C) ];

It is a simple circuit consisting of a capacitor and a coil. The capacitor is charged to voltage U0. The coil contains iron and has power dependent inductance; L = L0 / (1 + I^2). At time t = 0, the circuit is closed and the oscillation process is determined by two relationships: U = L(dI/dt) (voltage across the inductor) and I = -C(dU/dt) (current through the capacitor).

The differential equation at the beginning (du) is deduced from the two expressions above.

At t = 0, i = 0 and dI / dt = U0/L0. The solution I(t) is a periodic function which is more or less sinusoidal depending of how the U0 value is selected.

Given data: L0 = 1.10 H, C = 0.85μF. I've been trying four values of U0, from 300 V when the iron core influence is almost negligible, then 1200, 1800 and 2700 V where the power-curve is not sinusoidallike any longer.

I have assessed the magnitude of the oscillation-time T to

tfinal = 2*pi*sqrt(Lo*C);

for a circuit with constant C and constant L = L0 at T = 2π(sqrt(L0C)).

The task is to figure out a good algorithm to determine the peak of Imax with good accuracy, and a good algorithm to obtain the oscillation period T. I did it like this:

% Given start values

U      = 300;                      % Voltage
Lo     = 1.1;                      % Start value
C      = 0.85e-6;                  % Start value
tfinal = 2*pi*sqrt(Lo*C);          % Timespan

% ODE45

du = @(t,u)[ u(2); (2*u(1)/(1+u(1)^2))*u(2)^2-(1+u(1)^2)*u(1)/(Lo*C) ];       

opt = odeset('Reltol', 1e-10);                       
[tode,u] = ode45(du, [0 tfinal], [0 U/Lo], opt);    % Solving non-stiff ode
I = u(:,1);                                         

subplot(2,1,1), plot(tode,I,'b')                    % Plots vector with solutions
title([num2str(U) ' V']);
xlabel('Tid [s]');
ylabel('Ström [A]');
hold on
dI = u(:,2);                                        % Derivate of I

% Calculating the top value, Imax 

for b = 1:length(dI) - 1;                               
    if dI(b)*dI(b+1) >= 0 % Stepping to Imax due to derivate change from + to -
        b = b+1;                                      
    else
        disp([' Imax:'])
        disp([I(b)]);                               
        break
    end
end;

% Calculating the period T 

for n = 1:length(I) - 1;
    if I(n)*I(n+1) >= 0       % Stepping to where I changes from + to -
        n = n+1;
    else
        g = (I(n+1)-I(n))/(tode(n+1)-tode(n));
        T = 2*(tode(n)-(I(n)/g));    % T (x2 for the full period)
        disp([' Period time [seconds]:'])
        disp([T]);
        break
    end;
end;

The thing is that this works perfectly fine, but my teacher want's me to try to find the Imax and T values using RK4 and interpolation instead of my "stepping"-method. I'm seriously stuck here so any help would be truly and deeply appreciated!

What I've done so far with the RK4 (simply just the way the RK4 works):

for i=1:(length(x)-1)                              

      k_1 = du(x(i),y(i));
      k_2 = du(x(i)+0.5*h,y(i)+0.5*h*k_1);
      k_3 = du((x(i)+0.5*h),(y(i)+0.5*h*k_2));
      k_4 = du((x(i)+h),(y(i)+k_3*h));

      y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  
  end

Thank you in advance.