Ciao,
premetto che non ho il tempo per mettermi a controllare tutti i dettagli, ma ho dato un'occhiata.
Mi pare di vedere che tu sia interessato alla risoluzione di un problema del tipo $A \frac{dm}{dt}= F(\vec{m})$, dove $F$ è ovviamente una funzione a valori vettoriali. $A$ viene chiamata solitamente "matrice di massa", e ci sono solver dedicati a questo tipo di problemi in MatLab/Octave (tipo ode23t).
In particolare, l'inversione della matrice e il passaggio da $A(\vec{m}) \frac{dm}{dt}=F(m)$ a $\frac{dm}{dt}=A(m)^{-1} F(m)$, tramite il comando backslash, non è in genere una buona scelta. In questo caso inoltre la tua matrice dipende esplicitamente da $m=(m_1, \ldots,m_6)$.
Se non hai molto tempo, una strategia potrebbe essere quella di leggere la documentazione di un solver per problemi con matrici di massa, e passare la funzione e la matrice del problema in esame, che è una scemata ovviamente perché fai fare tutto a MatLab, ma avresti la sicurezza dell'implementazione.
Se vuoi farlo da solo e bene, allora la faccenda si complica perché un'implementazione naive di quel tipo non è detto sia corretta
Per la documentazione, leggi
qui, e questa utile pagina "choose an ODE solver":
clickBada che il metodo di Heun che tu proponi, è comunque un metodo esplicito, e quindi di facile implementazione
1. A tal proposito, mi sembra che il tuo codice possa essere migliorato. In particolare, io fossi in te, definirei la funzione $F$ una volta per tutti. Ossia
- Codice:
F=@(m) [m(1) - m(2);m(3)-m(5);m(6)-m(2);m(3);m(1)-m(4);]
dove ovviamente ho messo valori a caso in ogni entrata. La parte noiosa è quella di inserire in ognuna della sei componenti i vari prodotti vettoriali che dipenderanno da $m$, ma con comandi analoghi a quelli fatti da te. Insomma, alla fine ti ritroverai con una funzione vettoriale $F(\vec{m})$ ben precisa, un termine noto $m_0$ che definirai come un vettore a
sei componenti che è noto e ti ritroverai a risolvere il problema
$A(\vec{m}) \frac{dm}{dt}=F(\vec{m})$
con dato iniziale $\vec{m_0}=[m_1,\ldots,m_6]$.
Bada bene che $A$ è una matrice "state dependent", quindi per l'eventuale inversione (a.k.a risoluzione di un sistema lineare ad ogni time-step) dovresti sapere a priori che questa è non singolare
ad ogni time step, e questo può essere un problema. Ad ogni modo, seguendo il tuo ragionamento, ti riconduci al problema $\vec{m}'=A^{-1}(\vec{m}) F(\vec m)= \tilde{F}(\vec{m})$ e applichi Heun al problema $\vec{m}'=\tilde{F}(\vec{m})$, cioè
$y_{n+1}=y_n + \frac{h}{2}(\tilde{F}(y_n) + \tilde{F}(y_n+h \tilde{F}(y_n)))$
che con un banalissimo ciclo for si risolve.
Tirando lo somme, io farei così:
1. Prima, se hai effettivamente bisogna dei risultati, usa un solutore di MatLab. Sicuramente la tua soluzione sarà meno accurata, senza contare che MatLab utilizza un passo variabile,e mille altre cose che rendono il codice realemente efficiente.
2. Prova la risoluzione che ti ho proposto, che altro non è che la tua idea scritta in modo più schematizzato. A mio avviso, il tuo codice è un po' difficile da leggere, nel senso che quel ciclo per il metodo di Heun, facendo come l'ho proposto io, risulta effettivamente chiaro, e tutta la complessità della scrittura del codice viene mappata nello scrivere la funzione $F(\vec{m})$ una volta per tutte, ma non è la fine del mondo.
In bocca al lupo