Runge-Kutta 4 implemetation blowing up

Hello.

I am trying to solve the ODE system shown bellow using Runge-Kutta melhod, to compare with the solution given by the function ode23tb, that I successfully used to solve the exact same problem. I double checked the calculations and they seem to be correct. How can I avoid the extremely large numbers that I am getting in the Runge-Kutta coeficients?

fc = 10e3;
h = 1/(2000*fc);
t = h:h:0.01;

vD = 3;
vG = 10;
iL = 10/6.6;

s = double(...
    [vD*ones(1,length(t));
     vG*ones(1,length(t));
     iL*ones(1,length(t))]);
 
for i = 1:length(t)-1
    k1 = f(t(i),s(:,i));
    k2 = f(t(i)+h/2,s(:,i)+h*k1/2);
    k3 = f(t(i)+h/2,s(:,i)+h*k2/2);
    k4 = f(t(i)+h,s(:,i)+h*k3);

    s(:,i+1) = s(:,i+1) + 1/6*(k1 + 2*k2 + 2*k3 + k4);
end
plot(s(1,:));


function dydt = f(t,y)
    VDD = 10;
    Ld = 120e-6;
    Cg = 50e-12;
    Cb = 40e-12;
    RS = 50;
    RL = 6.6;
    
    vG = y(1); vD = y(2);
    iD = 10*(1/2.*(vG-3)+1/20.*log(2*cosh(10*(vG-3)))).*(1+0.003.*vD).*tanh(vD);
    vs = 3+20*sin(2*pi*10e3*t);
    dydt = [
            1/Cg*((vs-y(1))/RS-iD-y(3)-y(2)/RL);
            1/Cg*(vs-y(1))/RS-(Cb+Cg)/(Cb*Cg)*(iD+y(3)+y(2)/RL);
            (y(2)-VDD)/Ld
           ];      
end

[David Wilson]

It seems you have quite a stiff system to solve, and you are solving it for a very, very long time. You have two options:

(1) Use a decent stiff integrator such as ode23s, and let it choose the step size. Even then, you will need a small step size, and it seems no real point to go further than tfinal = 3e-4. The rest is just repetition.

(2) If you insist on using the brain-dead RK4, then you will need an extremely small step size to prevent instability. For example ode45 required h in the order of 1e-10. Of course then, you have no realistic idea of the accuracy, even if it is stable.

[José Filipe]

That seems indeed to be the problem. Could you please write your comment as an answer so that I can accept your answer?

Answer 1

[James Tursa]

s(:,i+1) = s(:,i) + 1/6*(k1 + 2*k2 + 2*k3 + k4); % changed s(:,i+1) to s(:,i) on rhs

[José Filipe]

Hi! Thank you very much for your answer and for finding an error. However, this seems not to be the sole problem. It stil fails to compute the first k2, since it goes to infinity. (As a side note, I fixes another error that i did while traying to scale down the coeficients).