Demande d'aide pour programmation matlab

Bonjour,

Je suis étudiante en licence de physique , et débutante avec Matlab.

Le problème que j'ai à résoudre porte sur la résolution d'équation différentielle du couplage de 2 circuits RLC par condensateur sur Matlab par la méthode de Runge Kutta d'orde 4.

Voici à quoi ressemble le montage, l'énoncé et les équas diff de ce sujet : http://subaru2.univ-lemans.fr/enseignements/physique/02/electri/coupcond.html

Voici ensuite mon programme de résolution :

A=2e-7; (avec A=L*C1)
B=2e-6; ( B=R*C1)
C=1.04; ( C=1+C1/C0)
D=0.04; ( D=C1/C0)
E=0.14; (E=C1/C2 + C1/C0)
dt=0.1;
t=0:dt:10;
n=length(t);

%Conditions initiales

q1=1e-5; (q1 correspond à la charge de C1)
q2=0; (q2 la charge de C2)
q3=50e-3; (q3 correspond à dq1/dt,donc à I1)
q4=0; (q4 à dq2/dt, c'est I2)

q(:,1)=[q1 q2 q3 q4]';
matrice=[0 0 1 0 ;
0 0 0 1 ;
(-C/A) D (-B/A) 0 ;
(D/A) (-E/A) 0 (-B/A)]

%Résolution par RK4

for i=1:n-1
K1=dt*fonc(q(:,i),matrice);
K2=dt*fonc(q(:,i)+K1/2,matrice);
K3=dt*fonc(q(:,i)+K2/2,matrice);
K4=dt*fonc(q(:,i)+K3,matrice);

q(:,i+1)=q(:,i)+(1/6)*(K1+2*K2+2*K3+K4);

t(i+1)=i*dt;
end

plot(t,q(1:2,:)) (Je cherche a tracer I1 et I2 en fonction de t)



Et pour finir, voici ma fonction:

function f=fonc(q,matrice);
f=matrice*q



Le problème est que ce programme ne marche pas. J'obtiens des valeurs hallucinantes dans la boucle, qui vont jusqu'a du e+302 pour les composantes de ma fonction, qui finit par m'afficher Not A Number à la fin! Sans parler du graphe qui ne ressemble pas....à celui d'un circuit RLC!


Si quelqu'un pouvait m'aider...Très franchement, ça ne serait pas de refus, car pour le moment je patauge lamentablement!