\documentclass{template-jurnal} %Bagian ini diedit oleh editor Jurnal Matematika UNAND
\usepackage{float}
\hyphenation{di-tulis-kan di-terang-kan}
%Perhatikan aturan penulisan dan ukuran huruf yang digunakan
\begin{document}

\markboth{Susila Bahri dkk} %Jika lebih dari dua penulis, tuliskan sebagai Nama Penulis Pertama dkk.
{Analisis Kestabilan Model Dinamika Perceraian \textit{MVQEDR}}

%%%%%%%%%%%%%%%%%%%%% Publisher's Area please ignore %%%%%%%%%%%%%%%
%
\catchline{}{}{}{}{}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{ANALISIS KESTABILAN MODEL DINAMIKA PERCERAIAN \textit{MVQEDR}}

\author{SUSILA BAHRI$^{a}$\footnote{penulis korespondensi}, MIYA QARLINA HUTAGALUNG$^{b}$, EFENDI$^{c}$, MUHAFZAN$^{d}$\\ }

\address{$^{a}$ Universitas Andalas,\\
$^{b}$ Universitas Andalas,\\
$^{c}$ Universitas Andalas.\\
$^{d}$ Universitas Andalas.\\
email : \email{susilabahri@sci.unand.ac.id, 2010431012\_miya@student.unand.ac.id, efendi@sci.unand.ac.id, muhafzan@sci.unand.ac.id}}

\maketitle
\setcounter{page}{1} %bagian ini diedit oleh editor Jurnal Matematika UNAND

\begin{abstract}
\begin{center}
Diterima ..... \quad Direvisi ..... \quad Dipublikasikan ..... %tanggal-tanggal tersebut \textbf{dikosongkan} saja
\end{center}

\bigskip

\textbf{Abstrak}. %Dalam bahasa Indonesia
Tiga faktor penyebab perceraian, yaitu masalah ekonomi rumah tangga, perselisihan dan pertengkaran terus- menerus, dan kekerasan dalam rumah tangga, masih memberikan kontribusi besar terhadap angka perceraian di Indonesia. Meskipun peme-\breakline rintah telah melakukan upaya penyuluhan bagi ketiga kelompok rumah tangga tersebut, namun pada kenyataannya kasus perceraian tidak kunjung berkurang. Oleh karena itu, perlu diketahui secara pasti seberapa besar pengaruh penyuluhan yang telah dilaksanakan oleh pemerintah terhadap kelompok ini. Pada penelitian ini, model matematika $MVQEDR$ terlebih dahulu dibangun. Analisis kestabilan titik ekuilibrium model dilakukan dengan menentukan nilai eigen dan matriks Jacobian dan diperoleh bahwa titik ekuilibrium bebas perceraian  stabil asimtotik jika $R_0 =  0,003111368984 < 1$ dan titik ekuilibrium endemik tidak stabil asimtotik jika $R_0 =  1,065035325 > 1$. Simulasi numerik dilakukan dengan menggunakan perangkat lunak MAPLE.

\bigskip

\textbf{Abstract}. % Dalam bahasa Inggris
\textit{The three factors that cause divorce, namely household economic pro-\breakline blems, continuous disputes and quarrels, and domestic violence, still contribute greatly to the divorce rate in Indonesia. Although the government has conducted counseling efforts for these three groups of households, in reality divorce cases have not decreased. Therefore, it is necessary to know exactly how much influence the counseling that has been implemented by the government has on this group. In this study, a mathematical model of $MVQEDR$ is first constructed. The stability analysis of the equilibrium point of the model is carried out by determining the eigenvalues and Jacobian matrix and it is obtained that the divorce-free equilibrium point is asymptotically stable if $R_0 = 0,003111368984 < 1$ and the endemic equilibrium point is not asymptotically stable if $R_0 = 1,065035325 > 1$. Numerical simulations were carried out using MAPLE software.}

\end{abstract}

\keywords{Cerai, Model Matematika, Simulasi Numerik.}

\section{Pendahuluan}

Perceraian adalah putusnya atau pembubaran perkawinan antara dua orang \cite{A}. Berdasarkan data Badan Pusat Statistik (BPS) tahun 2024, meskipun jumlah kasus perceraian di Indonesia berfluktuasi, namun jumlah kasus tersebut cenderung meningkat. Peningkatan ini disebabkan oleh beberapa faktor yaitu poligami, cacat tubuh, kekerasan dalam rumah tangga (KDRT), perselisihan dan pertengkaran terus-menerus, masalah ekonomi, dan zina. Di antara faktor tersebut, masalah ekonomi, perselisihan dan pertengkaran terus-menerus, dan Kekerasan Dalam Rumah Tangga (KDRT) selalu mendominasi dalam 5 tahun terakhir \cite{B}.
De-\breakline ngan kata lain, tiga kelompok rumah tangga dengan faktor tersebut telah menjadi penyumbang terbesar peningkatan jumlah kasus perceraian di Indonesia.

Banyak usaha atau program yang telah dilakukan pemerintah demi menekan angka perceraian antara lain membatasi perkawinan dini, kampanye tentang pen-\breakline tingnya keharmonisan rumah tangga, dan kewajiban mengikuti sesi konseling (\textit{counseling}) saat proses perceraian. Bagaimanapun juga diantara usaha tersebut,  konse-\breakline ling menjadi hal yang sangat penting dan menentukan karena konseling berhubu-\breakline ngan dengan usaha penanganan atau penyelesaian langsung terhadap masalah rumah tangga pada saat proses perceraian \cite{C}. Meskipun demikian, usaha konseling yang telah dilaksanakan pemerintah terlihat belum dapat menekan angka perceraian. Karenanya, perlu dikaji seberapa besar pengaruh konseling dan intensitas konseling yang telah dilakukan pemerintah selama ini, khususnya terhadap tiga kelompok rumah tangga penyumbang terbesar kasus perceraian tersebut. Hal ini penting diketahui agar konseling dapat menjadi program yang lebih efektif dalam mencegah perceraian sehingga jumlah kasus perceraian dapat diminimumkan. Pada penelitian ini juga dikaji, seberapa besar pengaruh individu rumah tangga yang tidak melakukan konseling terhadap peningkatan jumlah kasus perceraian di Indonesia.

Penelitian tentang perceraian yang menganalogikan perceraian sebagai penyakit menular telah dilakukan oleh beberapa peneliti. Patience Pokuaa Gambrah dkk (2018) menguji pengaruh konseling terhadap kesehatan kelompok individu yang mengalami perceraian. Penelitiannya menunjukkan bahwa konseling dapat meningkatkan kesehatan subpopulasi penderita perceraian \cite{D}. Reubeu Iorfyer Gwe-\breakline ryina dkk (2021) mengkonstruksi model epidemi perceraian dengan menggunakan fungsi insidensi standar. Hasil numerik menyatakan bahwa mediasi dan rekonsiliasi dapat menstabilkan perkawinan \cite{A}. Haileyesus Tessema dkk (2022) mendapatkan penularan perceraian disebarkan oleh perempuan atau laki-laki yang bercerai setelah menikah. Peneliti berhasil menentukan nilai parameter  penularan perceraian di Etiopia \cite{E}. Berbeda dengan yang dilakukan oleh para peneliti tersebut, pada tulisan ini dibahas tentang peningkatan perceraian ditinjau dari segi faktor-faktor penyebab terbesar perceraian di Indonesia, khususnya KDRT, perselisihan dan pertengkaran terus-menerus, dan masalah ekonomi.

\section{Metode}

Pertama, model matematika dinamika perceraian $MVQEDR$ dibangun dengan mengelompokkan populasi ke dalam enam kompartemen yaitu \textit{\textbf{M}} adalah kelompok individu yang sudah menikah, \textit{\textbf{V}} adalah kelompok individu yang mengalami KDRT, \textit{\textbf{Q}} adalah kelompok individu yang mengalami perselisihan dan pertengkaran terus-menerus, \textit{\textbf{E}} adalah kelompok individu yang mengalami masalah ekonomi rumah tangga, \textit{\textbf{D}} adalah kelompok individu yang bercerai, dan \textit{\textbf{R}} adalah kelompok individu yang sembuh/sadar. Parameter $\pi, \alpha, \beta, \gamma, \delta, C, \epsilon, \zeta, \theta, \sigma, \rho, \psi, \omega$
secara berturut-turut merupakan laju kelahiran alami, laju kematian alami, laju individu menikah kontak dengan individu bercerai dan berpindah ke kelompok yang meng-\breakline alami KDRT, laju individu menikah kontak dengan individu bercerai dan berpindah ke kelompok yang mengalami perselisihan dan pertengkaran terus-menerus, laju individu menikah kontak dengan individu bercerai dan berpindah ke kelompok yang mengalami masalah ekonomi  rumah tangga, laju individu yang melakukan konse-\breakline ling, laju perpindahan individu yang mengalami KDRT ke kelompok sembuh/sadar, laju perpindahan individu yang mengalami KDRT ke kelompok bercerai, laju perpindahan individu yang mengalami perselisihan dan pertengkaran terus-menerus ke kelompok sembuh/sadar, laju perpindahan individu yang mengalami pertengkaran dan perselisihan terus-menerus ke kelompok bercerai, laju perpindahan individu yang mengalami masalah ekonomi  rumah tangga ke kelompok sembuh/sadar, laju perpindahan individu yang mengalami masalah ekonomi  rumah tangga ke kelompok bercerai, dan laju perpindahan individu sembuh/sadar ke kelompok menikah. 

Hubungan antar kompartemen digambarkan pada Gambar 1.
Titik ekuilibrium bebas perceraian ditentukan dengan membuat subpopulasi bercerai sama dengan nol $(D=0)$. Titik ekuilibrium endemik diperoleh dengan mengasumsikan bahwa subpopulasi bercerai $(D>0)$. Kemudian untuk menentukan bilangan reproduksi dasar digunakan metode Next Generation Matrix \cite{F}. Selanjutnya, jenis kestabilan dari kedua titik ekuilibrium ditentukan dengan menggunakan nilai eigen dan matriks Jacobian \cite{G}. Pada tahap akhir, dilakukan simulasi numerik yang menggambarkan dinamika masing-masing populasi dengan menggunakan \textit{software} MAPLE.

\section{Hasil dan Pembahasan}
\subsection{Model Matematika}
Populasi model dibagi menjadi enam kategori: individu menikah (\textit{Susceptible}), di-\breakline lambangkan dengan $M$; individu yang mengalami kekerasan dalam rumah tangga (\textit{Exposed}), dilambangkan dengan $V$, individu yang mengalami perselisihan terus-menerus (\textit{Exposed}), dilambangkan dengan $Q$, individu yang mengalami masalah ekonomi rumah tangga (\textit{Exposed}), dilambangkan dengan $E$, individu bercerai (\textit{Infected}), dilambangkan dengan $D$, dan individu sembuh/sadar (\textit{Recovered}), dilambangkan dengan $R$. Individu dari kelompok rentan dapat berpindah ke kelas terpapar jika mereka berinteraksi dengan kelas terinfeksi. Selanjutnya, kelas terpapar dapat berpindah ke kelas terinfeksi jika mereka tidak mengikuti konseling. Kelas terpapar dapat pindah ke kelas sembuh/sadar jika mereka melakukan konseling. Kelas sembuh/sadar dapat berpindah ke kelas rentan jika mereka kembali ke pernikahan mereka. Populasi dilambangkan dengan $N$. 
Berdasarkan kasus tersebut, diagram model matematika perceraian dapat digambarkan sebagai berikut: 
\clearpage
\begin{figure}[htbp]
\center{\includegraphics[width=8cm]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/BaganMVQEDR.jpg}}
 \caption{Diagram Model $MVQEDR$} \label{gbr1}
\end{figure}

Diagram pada Gambar 1 dapat dikonstruksi ke dalam sistem persamaan dife-\breakline rensial nonlinier berikut ini:
\begin{equation}
\resizebox{!}{2.cm}{$
\begin{aligned}
\frac{dM}{dt} &= \pi+\omega R-M\left(\alpha+\beta D+\gamma D+\delta D\right),\\
\frac{dV}{dt} &= \beta MD-V\left(\alpha+\epsilon C+\zeta (1-C)\right),\\
\frac{dQ}{dt} &= \gamma MD-Q\left(\alpha+\theta C+\sigma (1-C)\right),\\
\frac{dE}{dt} &= \delta MD-E\left(\alpha+\rho C+\psi(1-C)\right),\\
\frac{dD}{dt} &= \zeta(1-C) V +\sigma(1-C) Q+\psi(1-C)E-\alpha D,\\
\frac{dR}{dt} &= \epsilon CV + \theta CQ + \rho CE-R\left(\alpha+\omega\right).
\end{aligned}$} \label{kompartemen}
\end{equation}
dengan kondisi awal pada $t = 0$ adalah:
$M(0) \geq 0, \quad V(0) \geq 0, \quad Q(0) \geq 0, \quad E(0) \geq 0, \quad D(0) \geq 0, \quad R(0) \geq 0$. Dalam hal ini, parameter-parameter $\pi, \alpha, \beta, \gamma, \delta, C, \epsilon, \zeta, \theta, \sigma, \rho, \psi, \omega$ adalah konstanta positif. Deskripsi dari masing-masing parameter dapat dilihat pada Tabel I berikut.
\vspace{-0.5cm}
\begin{table}[h!]
  \centering
  \caption{Nilai Parameter Subpopulasi Perceraian Tahun 2019-2023}
  \label{tbllvariabel}
  \resizebox{!}{1.7cm}{$\begin{tabular}{|c|c|c||c|c|c|}
    \hline
    \textbf{Parameter} & \textbf{Nilai}  & \textbf{Referensi} & \textbf{Parameter} & \textbf{Nilai} & \textbf{Referensi} \\
    \hline
    $\alpha$ & 0,00166& \cite{B} &$\theta$ & 0,1 & Asumsi\\
    \hline
    $\pi$ & 0,0033 & \cite{B}  & $\rho$ & 0,1 & Asumsi\\
    \hline
    $\beta$ & 0,000592 & \cite{B} & $\zeta$ & 0,000001 & Asumsi\\
    \hline
    $\gamma$ & 0,0278 & \cite{B}  & $\sigma$ & 0,000001 & Asumsi\\
    \hline
    $\delta$ & 0,01194 & \cite{B}  & $\psi$ & 0,000001 & Asumsi\\ \hline
    $C$ & 0,12 & Asumsi & $\omega$ & 0,2 & Asumsi\\ \hline
    $\epsilon$ &  0,1 & Asumsi & & &  \\ \hline
    \end{tabular}$} \label{tabelparameterbbs}
\end{table}
\subsection{Analisis Model}
\subsubsection{Titik Ekuilibrium}
Titik-titik ekuilibrium dari sistem (\ref{kompartemen}) dapat diperoleh dengan mengatur setiap persamaan menjadi nol \cite{H} \cite{I}:
\begin{equation}
\resizebox{!}{1.5cm}{$
\begin{aligned}
0 &= \pi+\omega R-M\left(\alpha+\beta D+\gamma D+\delta D\right)  \\ 
0 &= \beta MD-V\left(\alpha+\epsilon C+\zeta (1-C)\right)    \\
0 &= \gamma MD-Q\left(\alpha+\theta C+\sigma (1-C)\right) \\
0 &= \delta MD-E\left(\alpha+\rho C+\psi(1-C)\right)  \\
0 &= \zeta(1-C) V +\sigma(1-C) Q+\psi(1-C)E-\alpha D \\
0 &= \epsilon CV + \theta CQ + \rho CE-R\left(\alpha+\omega\right) .
\end{aligned}$}\label{MVQEDR} 
\end{equation}
Dua titik ekuilibrium yaitu titik ekuilibrium bebas perceraian dan titik ekuilibrium endemik perceraian untuk sistem (\ref{MVQEDR}) adalah
\begin{enumerate}
    \item Titik Ekuilibrium Bebas Perceraian\\
Titik ekuilibrium ini diperoleh dengan mengasumsikan $D=0$ sehingga $E_T^0=(M^0,V^0,Q^0,E^0,D^0,R^0)= \left(\frac{\pi}{\alpha}, 0, 0, 0, 0, 0\right)$.
\item Titik Ekuilibrium Endemik\\
Titik ekuilibrium endemik yang mengasumsikan adanya kasus perceraian $D>0$ adalah  $E_T^*=(M^*,V^*,Q^*,E^*,D^*,R^*)$ dengan
\vspace{-0.01cm}
\begin{equation*}
\resizebox{!}{3cm}{$
\begin{aligned}
M^* &= \frac{a}{b} \notag \\
V^* &= k_1 \left(\frac{\left(\frac{\pi}{\alpha}-\frac{a}{b}\right)}{1+k_1+k_2+k_3+k_4} \right) \notag \\
Q^* &= k_2  \left(\frac{\left(\frac{\pi}{\alpha}-\frac{a}{b}\right)}{1+k_1+k_2+k_3+k_4}\right) \notag \\
E^* &= k_3  \left(\frac{\left(\frac{\pi}{\alpha}-\frac{a}{b}\right)}{1+k_1+k_2+k_3+k_4}\right)  \label{E*} \\
D^* &= \frac{\left(\frac{\pi}{\alpha}-\frac{a}{b}\right)}{1+k_1+k_2+k_3+k_4}  \notag \\
R^* &= k_4  \left(\frac{\left(\frac{\pi}{\alpha}-\frac{a}{b}\right)}{1+k_1+k_2+k_3+k_4}\right), 
\end{aligned}$}
\end{equation*}
dan
\begin{equation*}
\resizebox{!}{2.5cm}{$
\begin{aligned}
a&= \alpha \\
b&= \frac{\zeta(1-C) \beta}{\alpha + \epsilon C + \zeta (1-C)} + \frac{\sigma(1-C) \gamma}{\alpha + \theta C + \sigma (1-C)} + \frac{\psi(1-C) \delta}{\alpha + \rho C + \psi(1-C)}\\
 k_1&= \frac{a \beta }{b(\alpha + \epsilon C + \zeta (1-C))}   \\ 
k_2&= \frac{a \gamma }{b(\alpha + \theta C + \sigma (1-C))} \\
k_3&= \frac{a \delta }{b(\alpha + \rho C + \psi(1-C))}  \\
k_4&= \frac{\frac{a \epsilon C \beta  }{b(\alpha + \epsilon C + \zeta (1-C))} +   \frac{a \theta C \gamma}{b(\alpha + \theta C + \sigma (1-C))} +  \frac{a \rho C \delta }{b(\alpha + \rho C + \psi(1-C))}}{(\alpha + \omega)}.
\end{aligned}$}
\end{equation*}

\subsubsection{Bilangan Reproduksi Dasar ($R_0$)}
Bilangan Reproduksi Dasar dapat ditentukan dengan mengaplikasikan metode \textit{Next Generation Matrix} pada persamaan kompartemen terinfeksi berikut \cite{J} :
\begin{equation}
\resizebox{!}{1.5cm}{$
\begin{aligned}
    \frac{dV}{dt} &= \beta MD-V\left(\alpha+\epsilon C+\zeta (1-C)\right),\\
    \frac{dQ}{dt} &= \gamma MD-Q\left(\alpha+\theta C+\sigma (1-C)\right),\\
    \frac{dE}{dt} &= \delta MD-E\left(\alpha+\rho C+\psi(1-C)\right),\\
    \frac{dD}{dt} &= \zeta(1-C) V +\sigma(1-C) Q+\psi(1-C)E-\alpha D
\end{aligned}$} \label{terinfeksi}
\end{equation}
Langkah pertama adalah menentukan matriks transmisi F dan matriks transisi V, dan radius spektral $K=FV^{-1}$ sehingga diperoleh \cite{L}
\begin{equation*}
\resizebox{!}{1cm}{$
\begin{aligned}
   F &= \begin{bmatrix}
       0 & 0 & 0 & \frac{\beta \pi}{\alpha}\\
       0 & 0 & 0 & \frac{\gamma \pi}{\alpha}\\
       0 & 0 & 0 & \frac{\delta \pi}{\alpha}\\
       0 & 0 & 0 & 0
   \end{bmatrix}.
   \end{aligned}$}
\end{equation*}   
Kemudian,
\begin{equation*}
\resizebox{!}{0.8cm}{$ 
\begin{aligned}
   V &= \begin{bmatrix}
       \alpha+\epsilon C+\zeta (1-C)  & 0 & 0 & 0\\
       0 & \alpha+\theta C+\sigma (1-C) & 0 & 0\\
       0 & 0 &\alpha+\rho C+\psi(1-C)  & 0\\
        -\zeta(1-C)  & -\sigma(1-C) & -\psi(1-C) & \alpha 
   \end{bmatrix}.
\end{aligned}$}
\end{equation*}
Sehingga, \textit{Next Generation Matrix}  $K=FV^{-1}$ adalah
\begin{align*}
K &= FV^{-1}\\
&= \begin{bmatrix}
0 & 0 & 0 & \frac{\beta \pi}{\alpha}\\
0 & 0 & 0 & \frac{\gamma \pi}{\alpha}\\
0 & 0 & 0 & \frac{\delta \pi}{\alpha}\\
0 & 0 & 0 & 0
\end{bmatrix} \begin{bmatrix}
\frac{1}{ \alpha+\epsilon C+\zeta (1-C)} & 0 & 0 & 0 \\
0 & \frac{1}{\alpha+\theta C+\sigma (1-C)} & 0 & 0 \\
0 & 0 & \frac{-1}{\alpha+\rho C+\psi(1-C)} & 0 \\
\frac{\zeta(1 - C)}{ \alpha(\alpha+\epsilon C+\zeta (1-C))} & \frac{\sigma(1 - C)}{\alpha(\alpha+\theta C+\sigma (1-C))} & \frac{\psi(1 - C)}{\alpha(\alpha+\rho C+\psi(1-C))} & \frac{1}{\alpha}
\end{bmatrix}
\end{align*}
\begin{align*}
&= \begin{bmatrix} \frac{\beta \pi \zeta(1-C)}{\alpha^2(\alpha + \epsilon C + \zeta(1-C))} & \frac{\beta \pi \sigma(1-C)}{\alpha^2(\alpha + \theta C + \sigma(1-C))} & \frac{\beta \pi \psi(1-C)}{\alpha^2(\alpha + \rho C + \psi(1-C))} & \frac{\beta \pi}{\alpha^2} \\
\frac{\gamma \pi \zeta(1-C)}{\alpha^2(\alpha + \epsilon C + \zeta(1-C))} & \frac{\gamma \pi \sigma(1-C)}{\alpha^2(\alpha + \theta C + \sigma(1-C))} & \frac{\gamma \pi \psi(1-C)}{\alpha^2(\alpha + \rho C + \psi(1-C))} & \frac{\gamma \pi}{\alpha^2} \\
\frac{\delta \pi \zeta(1-C)}{\alpha^2(\alpha + \epsilon C + \zeta(1-C))} & \frac{\delta \pi \sigma(1-C)}{\alpha^2(\alpha + \theta C + \sigma(1-C))} & \frac{\delta \pi \psi(1-C)}{\alpha^2(\alpha + \rho C + \psi(1-C))} & \frac{\delta \pi}{\alpha^2} \\
0 & 0 & 0 & 0 \end{bmatrix}.
\end{align*}
Bilangan reproduksi dasar ($R_0$) adalah radius spektral atau nilai absolut maksimum dari nilai eigen \textit{Next Generation Matrix} K, sedemikian rupa sehingga:
\begin{align*} 
det\left(K - \lambda I\right)&=&0  \notag \\
det\begin{pmatrix} \frac{\beta \pi \zeta(1-C)}{\alpha^2(\alpha + \epsilon C + \zeta(1-C))}-\lambda & \frac{\beta \pi \sigma(1-C)}{\alpha^2(\alpha + \theta C + \sigma(1-C))} & \frac{\beta \pi \psi(1-C)}{\alpha^2(\alpha + \rho C + \psi(1-C))} & \frac{\beta \pi}{\alpha^2} \\
\frac{\gamma \pi \zeta(1-C)}{\alpha^2(\alpha + \epsilon C + \zeta(1-C))} & \frac{\gamma \pi \sigma(1-C)}{\alpha^2(\alpha + \theta C + \sigma(1-C))} -\lambda & \frac{\gamma \pi \psi(1-C)}{\alpha^2(\alpha + \rho C + \psi(1-C))} & \frac{\gamma \pi}{\alpha^2} \\
\frac{\delta \pi \zeta(1-C)}{\alpha^2(\alpha + \epsilon C + \zeta(1-C))} & \frac{\delta \pi \sigma(1-C)}{\alpha^2(\alpha + \theta C + \sigma(1-C))} & \frac{\delta \pi \psi(1-C)}{\alpha^2(\alpha + \rho C + \psi(1-C))} -\lambda & \frac{\delta \pi}{\alpha^2} \\
0 & 0 & 0 & -\lambda \end{pmatrix} &=&0 
\end{align*}
    $\lambda_1 = \frac{A+B+C+D}{E+F+G}$ dan
    $\lambda_{2,3,4} =0$, dengan 

\noindent
\begin{equation*}
\resizebox{!}{2cm}{$
\begin{aligned}
A &= C^3 \left( \zeta \beta (\psi \theta - \rho \theta) + \delta (\psi \theta - \psi \sigma) + \gamma (\rho \sigma - \psi \sigma) + \delta \epsilon (\psi \sigma - \psi \theta) + \epsilon \gamma (\psi \sigma - \rho \sigma) \right) \\
B &= C^2 ( \zeta \alpha (\beta (\sigma - \theta) + \delta \psi + \gamma \sigma) + \beta (3 \psi \sigma - 2 \psi \theta - 2 \rho \sigma + \rho \theta) + \delta (3 \psi \sigma - 2 \psi \theta) \\
& + \gamma (3 \psi \sigma - 2 \rho \sigma) - \alpha (\delta \epsilon (\psi - \sigma) + \epsilon \gamma \sigma - 2 \delta \psi \sigma + \gamma (\psi \sigma - \rho \sigma)) ) \\
C &= C \left( -\zeta \alpha^2 \beta - 2 \zeta \alpha (\beta \psi + \delta \psi + \gamma \sigma) - 3 \zeta (\psi \sigma + \delta \psi + \gamma \psi \sigma) + \alpha (\delta \epsilon \psi - \epsilon (\psi \sigma - \rho \sigma)) \right) \\
D &= \zeta (\alpha^2 \beta + \alpha (\beta \psi + \sigma) + \delta \psi + \gamma \sigma)+\alpha (\delta \psi \sigma + \gamma \psi \sigma) \\
E &= C^3 \left( \zeta (\psi \theta - \psi \sigma) + \epsilon (\psi \sigma - \rho \sigma) \right) + C^2 ( \zeta (\alpha (\sigma - \theta) + 3 \psi \sigma - 2 \psi \theta - 2 \rho \sigma + \rho \theta) \\
& - \alpha (\epsilon (\sigma - \theta) + 2 \psi \sigma - \psi \theta - \rho \sigma)) \\
F &= C \left( -\zeta \alpha^2 - 2 \zeta (\alpha \psi + \beta \sigma) - 3 (\psi \sigma + \delta \psi \sigma + \gamma \psi \sigma) + \alpha (\epsilon \psi - \psi \sigma + \rho \sigma) \right) \\
G &= \zeta (\alpha^2 + \alpha (\psi + \sigma) + \psi \sigma) + \alpha (\psi \sigma + \sigma). 
\end{aligned}$}  
\end{equation*}
Menurut \cite{I} \cite{J}, nilai eigen terbesar dari matriks $FV^{-1}$ adalah bilangan reproduksi dasar ($R_0$), sehingga
\begin{equation*}
\resizebox{!}{0.4cm}{$
\begin{aligned}
R_0 = \frac{A+B+C+D}{E+F+G}.
\end{aligned}$}
\end{equation*}
\subsection{Analisis Kestabilan Titik Ekuilibrium}
\subsubsection{Kestabilan Titik Ekuilibrium Bebas Perceraian}
Titik ekuilibrium bebas perceraian dari model (\ref{kompartemen}) stabil asimtotik jika nilai eigen matriks Jacobian negatif. Matriks Jacobian dari model (\ref{kompartemen}) pada titik ekuilibrium bebas perceraian ($J_{E^0})$) diberikan oleh
\begin{equation*}
\resizebox{!}{1.3cm}{$
\begin{aligned}
J_{E_T^0} = \begin{bmatrix}
g_1 & 0 & 0 & 0 & g_2 & g_3 \\
0 & g_4 & 0 & 0 & g_5 & 0 \\
0 & 0 & g_6 & 0 & g_7 & 0 \\
0 & 0 & 0 & g_8 & g_9 & 0 \\
0 & g_{10} & g_{11} & g_{12} & g_1 & 0 \\
0 & g_{13} & g_{14} & g_{15} & 0 & g_{16}
\end{bmatrix}.
\end{aligned}$}
\end{equation*}
Formula elemen matriks Jacobian titik ekuilibrium bebas perceraian disajikan dalam tabel berikut.

\begin{table}[h!]
\centering
\caption{Formula Elemen Matriks Jacobian Titik Ekuilibrium Bebas Perceraian}
 \label{formbebas}
\resizebox{!}{2.0cm}{$\begin{tabular}{|c|c|c|c|}
\hline
\textbf{Variabel} & \textbf{Formula} & \textbf{Variabel} & \textbf{Formula} \\ \hline
$g_1$  & $-\alpha$                                  & $g_9$   & $\frac{\delta \pi}{\alpha}$ \\ \hline
$g_2$  & $-\frac{\pi}{\alpha}(\beta + \gamma + \delta)$ & $g_{10}$  & $\zeta (1 - C)$ \\ \hline
$g_3$  & $\omega$                                   & $g_{11}$  & $\sigma (1 - C)$ \\ \hline
$g_4$  & $-(\alpha + \epsilon C + \zeta (1 - C))$   & $g_{12}$  & $\psi (1 - C)$ \\ \hline
$g_5$  & $\frac{\beta \pi}{\alpha}$                 & $g_{13}$  & $\epsilon C$ \\ \hline
$g_6$  & $-(\alpha + \theta C + \sigma (1 - C))$    & $g_{14}$  & $\theta C$ \\ \hline
$g_7$  & $\frac{\gamma \pi}{\alpha}$                & $g_{15}$  & $\rho C$ \\ \hline
$g_8$  & $-(\alpha + \rho C + \psi (1 - C))$        & $g_{16}$  & $-(\alpha + \omega)$ \\ \hline
\end{tabular}$} 
\end{table}
Persamaan karakteristik dapat diperoleh dengan menyelesaikan persamaan $\det(J_{E_T^0} - \lambda I)=0$ yaitu:

\begin{equation*}
\resizebox{!}{1.3cm}{$
\begin{aligned}
0&= \begin{vmatrix}
        g_1 - \lambda & 0 & 0 & 0 & g_2 & g_3 \\
    0 & g_4 - \lambda & 0 & 0 & g_5 & 0 \\
    0 & 0 & g_6 - \lambda & 0 & g_7 & 0 \\
    0 & 0 & 0 & g_8 - \lambda & g_9 & 0 \\
    0 & g_{10} & g_{11} & g_{12} & g_1 - \lambda & 0 \\
    0 & g_{13} & g_{14} & g_{15} & 0 & g_{16} - \lambda
    \end{vmatrix},
\end{aligned}$}
\end{equation*}
sehingga persamaan karakteristik yang sesuai adalah:
\begin{equation}
\resizebox{!}{0.3cm}{$
\begin{aligned}\label{uwu}
0&= (g_1 - \lambda)(g_{16} - \lambda)( \lambda^4 + a_1 \lambda^3 + a_2 \lambda^2 + a_3 \lambda + a_4),
\end{aligned}$}
\end{equation}
dengan
\begin{equation*}
\resizebox{!}{1.3cm}{$
\begin{aligned}
a_1 &= -(g_1 + g_4 + g_6 + g_8), \\
a_2 &= g_1 g_4 + g_1 g_6 + g_1 g_8 - g_{10} g_5 - g_{11} g_7 - g_{12} g_9 + g_4 g_6 + g_4 g_8 + g_6 g_8, \\
a_3 &= -g_1 g_4 g_6 - g_1 g_4 g_8 - g_1 g_6 g_8 + g_{10} g_5 g_6 + g_{10} g_5 g_8 + g_{11} g_4 g_7 + g_{11} g_7 g_8 \\
& + g_{12} g_4 g_9 + g_{12} g_6 g_9 - g_4 g_6 g_8, \\
a_4 &= g_1 g_4 g_6 g_8 - g_{10} g_5 g_6 g_8 - g_{11} g_4 g_7 g_8 - g_{12} g_4 g_6 g_9. 
\end{aligned}$}
\end{equation*}
Berdasarkan Definisi (2.2.3), titik ekuilibrium dikatakan stabil asimtotik jika $\lambda_i<0$.
Dari persamaan (\ref{uwu}), diperoleh 2 nilai eigen, yaitu
 $ \lambda_1 &= g_1 = -\alpha<0$ dan
 $\lambda_2 &= g_{16}= -(\alpha+\omega)<0$ .   

Dengan menggunakan kriteria Routh–Hurwitz, dapat dilihat bahwa bagian riil dari semua akar polinomial adalah negatif jika dan hanya jika $a_1 > 0, \frac{a_1 a_2}{a_3} > 1, \frac{a_3(a_1 a_2 - a_3)}{a_1^2 a_4} > 1$, dan $a_4 > 0$.

\subsubsection{Kestabilan Titik Ekuilibrium Endemik Perceraian}
Titik ekuilibrium endemik dari model (\ref{kompartemen}) stabil asimtotik jika nilai eigen matriks Jacobian negatif. Matriks jacobian model (\ref{kompartemen}) pada titik ekuilibrium endemik  $(J_{E^*})$ diberikan oleh
\begin{equation*}
\resizebox{!}{1.3cm}{$
\begin{aligned}
J_{E_T^*} = \begin{bmatrix}
h_1 & 0 & 0 & 0 & h_2 & h_3 \\
h_4 & h_5 & 0 & 0 & h_6 & 0 \\
h_7 & 0 & h_8 & 0 & h_9 & 0 \\
h_{10} & 0 & 0 & h_{11} & h_{12} & 0 \\
0 & h_{13} & h_{14} & h_{15} & h_{16} & 0 \\
0 & h_{17} & h_{18} & h_{19} & 0 & h_{20}
        \end{bmatrix}. \label{kestabilanE*}
\end{aligned}$}
\end{equation*}
Formula elemen matriks Jacobian titik ekuilibrium endemik perceraian disa-\breakline jikan dalam tabel berikut
\clearpage
\begin{table}[h!]
\centering
\caption{Formula Elemen Matriks Jacobian Titik Ekuilibrium Endemik Perceraian}
 \label{formendemik}
\resizebox{!}{2.8cm}{$\begin{tabular}{|c|c|c|c|}
\hline
\textbf{Variabel} & \textbf{Formula} & \textbf{Variabel} & \textbf{Formula} \\ \hline
$h_1$  & $-\alpha - \frac{(\beta + \gamma + \delta) \left(\frac{\pi}{\alpha} - \frac{a}{b}\right)}{1 + k_1 + k_2 + k_3 + k_4}$ 
& $h_{11}$ & $-(\alpha + \rho C + \psi (1-C))$ \\ \hline
$h_2$  & $-\frac{a(\beta + \gamma + \delta)}{b}$ 
& $h_{12}$ & $\frac{\delta a}{b}$ \\ \hline
$h_3$  & $\omega$ 
& $h_{13}$ & $\zeta (1-C)$ \\ \hline
$h_4$  & $\frac{\beta \left(\frac{\pi}{\alpha} - \frac{a}{b}\right)}{1 + k_1 + k_2 + k_3 + k_4}$ 
& $h_{14}$ & $\sigma (1-C)$ \\ \hline
$h_5$  & $-(\alpha + \epsilon C + \zeta (1-C))$ 
& $h_{15}$ & $\psi (1-C)$ \\ \hline
$h_6$  & $\frac{\beta a}{b}$ 
& $h_{16}$ & $-\alpha$ \\ \hline
$h_7$  & $\frac{\gamma \left(\frac{\pi}{\alpha} - \frac{a}{b}\right)}{1 + k_1 + k_2 + k_3 + k_4}$ 
& $h_{17}$ & $\epsilon C$ \\ \hline
$h_8$  & $-(\alpha + \theta C + \sigma (1-C))$ 
& $h_{18}$ & $\theta C$ \\ \hline
$h_9$  & $\frac{\gamma a}{b}$ 
& $h_{19}$ & $\rho C$ \\ \hline
$h_{10}$ & $\frac{\delta \left(\frac{\pi}{\alpha} - \frac{a}{b}\right)}{1 + k_1 + k_2 + k_3 + k_4}$ 
& $h_{20}$ & $-(\alpha + \omega)$ \\ \hline
\end{tabular}$} 
\end{table}

Persamaan karakteristik dapat diperoleh dengan menyelesaikan persamaan $\det(J_{E_T^*} - \lambda I)=0$ yaitu:
\begin{equation*}
\resizebox{!}{1.3cm}{$
\begin{aligned}
0 &=\begin{vmatrix}
h_1 - \lambda & 0 & 0 & 0 & h_2 & h_3 \\
h_4 & h_5 - \lambda & 0 & 0 & h_6 & 0 \\
h_7 & 0 & h_8 - \lambda & 0 & h_9 & 0 \\
h_{10} & 0 & 0 & h_{11} - \lambda & h_{12} & 0 \\
0 & h_{13} & h_{14} & h_{15} & h_{16} - \lambda & 0 \\
0 & h_{17} & h_{18} & h_{19} & 0 & h_{20} - \lambda
    \end{vmatrix}\notag \\
\end{aligned}$}
\end{equation*}
sehingga persamaan karakteristik untuk kestabilan titik ekuilibrium ini adalah 
\begin{align}\label{pkendemk}
0&= \lambda^6 + a_1 \lambda^5 + a_2\lambda^4+a_3 \lambda^3+a_4 \lambda^2+a_5 \lambda+a_6,
\end{align} 
dengan
\begin{equation*}
\scriptsize
\begin{aligned}
a_1 &= -h_1 - h_{16} - h_{20} - h_5 - h_{11} - h_8, \\
a_2 &= h_1(h_{11} + h_{16} + h_{20} + h_5 + h_8) + h_{11}(h_{16} + h_{20} + h_5 + h_8)- h_{12}h_{15} - h_{13}h_6 - h_{14}h_9 + h_{16}(h_{20} + h_5 + h_8)  \\
&\quad + h_{20}(h_5 + h_8) + h_5h_8, \\
a_3 &= -h_1(h_{16}h_{20} + h_{16}h_5 + h_{11}h_{16} + h_{20}h_5 + h_{11}h_{20} + h_{11}h_5)- h_{10}h_{15}h_2 - h_{16}(h_{20} + h_5)h_5 \\
&\quad  - h_{11}(h_{16}h_{20} + h_{16}h_5 + h_{20}h_5)+ h_{12}(h_{15}h_{20} + h_{15}h_5 + h_1h_{15}) - h_{10}h_{19}h_3 - h_1(h_{16} + h_{20} + h_5 + h_{11})h_8 \\
&\quad + h_{13}(h_6h_8 + h_1h_6 + h_{20}h_6 + h_{11}h_6) - h_{17}h_3h_4 - h_{13}h_2h_4, \\
a_4 &= h_{10}h_{16}h_{19}h_3 + h_{10}h_{19}h_3h_5 + h_1h_{16}(h_{20}h_5 + h_{5}h_8 + h_{11}h_{16})  + h_{11}(h_2h_{14}h_7 - h_{14}h_{20}h_9 - h_{14}h_5h_9)  \\
&\quad - h_{12}(h_{15}h_{20}h_8 + h_{15}h_5h_8) + h_{13}(h_2h_{20}h_4 + h_{11}h_{13}h_2h_4)+ h_{17}h_3h_4 - h_{11}h_{13}h_6 - h_{13}(h_{20} + h_6)h_6h_8 \\
&\quad - h_{11}(h_{13}h_6h_8 + h_{20}h_{13}h_6) + h_{10}h_{19}h_3h_8  + h_{1}h_{11}(h_{16}h_{20} + h_{16}h_5 + h_{20}h_5)  + h_{10}h_{15}(h_2h_{20} + h_2h_5) \\
&\quad - h_{12}(h_{15}h_{20} + h_{15}h_5), \\
a_5 &= h_{12}h_{15}h_{18}h_3h_7 - h_{14}h_{17}h_3h_6h_7 - h_{16}h_{18}h_3(h_5 + h_7) - h_1h_{11}(h_{16}h_{20} + h_{16}h_5 + h_{20}h_5)  \\
&\quad - h_{10}(h_{15}h_{2}h_{20}h_5 + h_{15}h_2h_{20}h_8)- h_{11}(h_{16}h_{20}h_5h_8 + h_{16}h_{20}h_8) + h_{12}(h_{15}h_{20}h_5 + h_{15}h_{20}h_8) \\
&\quad + h_{10}(h_{15}h_{18}h_3h_9 + h_{14}h_{19}h_3h_9) - h_{11}h_{16}h_{18}h_3h_7- h_{12}(h_{13}h_{19}h_3h_4 + h_{13}h_{18}h_3h_4) - h_{1}h_{11}h_{13}h_6h_8 \\
&\quad - h_{13}h_{20}h_6h_8 + h_{11}(h_{13}h_{20}h_6h_8 + h_{13}h_6h_8) - h_{11}h_{16}h_{17}h_3h_4 - h_{16}h_{17}h_3h_4h_8 - h_{11}h_{17}h_3h_4h_8, \\
a_6 &= -h_{12}h_{15}h_{18}h_3h_5h_7 + h_{10}h_{15}(h_{2}h_{20}h_5h_8 - h_{15}h_{17}h_3h_6)- h_{1}h_{12}h_{15}h_{20}h_5h_8 \\
&\quad  + h_{10}h_{15}(h_{2}h_{20}h_5h_8 + h_{15}h_{17}h_3h_6) - h_{11}(h_{13}h_{20}h_6h_8 + h_{13}h_{18}h_3h_6h_7) - h_{10}h_{15}(h_{17}h_{3}h_6 + h_{19}h_3h_5h_8)\\
&\quad  - h_{10}(h_{19}h_3h_8 - h_{16}h_{19}h_3h_8 + h_{11}h_{16}h_{20}h_5h_8).
\end{aligned}
\end{equation*}
Selanjutnya, menggunakan kriteria Routh–Hurwitz, dapat dilihat bahwa bagian riil dari semua akar polinomial adalah negatif jika dan hanya jika:

\begin{equation}\label{lop} 
\resizebox{!}{1.5cm}{$
\begin{aligned}
 & \quad a_1  > 0,  \\
 & \quad a_1 a_2 - a_3 > 0 , \\
& \quad a_3(a_1 a_2 - a_3) - a_1^2 a_4+a_1 a_5> 0 , \\ 
& \quad  a_1^2 (a_2 a_6 - a_4^2)  - a_1 r  + a_2 a_3 a_5 - a_3^2 a_4- a_5^2 > 0, \\
 & \quad -a_1^2 s - a_1 t  + a_3 u > 0,\\
& \quad v (w+x+y)>0, 
\end{aligned}$} 
\end{equation}
dengan 
\begin{equation*}\label{lopyy} 
\resizebox{!}{2.0cm}{$
\begin{aligned}
r&= a_2^2 a_5 - a_2 a_3 a_4 + a_3 a_6 - 2 a_4 a_5\\
s&=  a_1 a_6^2 - a_4^2 a_5 - a_3 a_4 a_6 + 2 a_2 a_5 a_6\\
t&= a_2^2 a_5^2 + a_2 a_3^2 a_6 - a_2 a_3 a_4 a_5 + 3 a_3 a_5 a_6 - 2 a_4 a_5^2\\
u&= a_2 a_5^2 + a_3^2 a_6 - a_3 a_4 a_5 - a_5^2 \\
v&= a_6\\
w&=- a_1^3 a_6^2+ a_1^2 ( 2 a_2 a_5 a_6 + a_3 a_4 a_6 - a_4^2 a_5)\\
x&=a_1 a_2 ( - a_2 a_5^2 + a_3 a_4 a_5 - a_3^2 a_6 ) + a_1 a_5 ( -3 a_3 a_6 + 2 a_4 a_5 )\\
y&=a_2 a_3 a_5^2 + a_3^2 ( a_3 a_6 - a_4 a_5 )- a_5^3.
\end{aligned}$}
\end{equation*}

\subsection{Simulasi Numerik}
Dengan menggunakan nilai untuk setiap parameter yang digunakan pada sistem (\ref{kompartemen}). Dari nilai parameter pada Tabel 1, diperoleh $R_0=0,003111368984< 1$ dan titik ekuilibrium bebas perceraian $E_T^0=(1,987951807; 0; 0; 0; 0, 0)$. Dalam simulasi titik ekuilibrium bebas perceraian, nilai awal digunakan $M(0)=0,0072594136, V(0)=0,000028809391, Q(0)=0,0008548875, E(0)=0,0004451261, D(0)=0,0016185539, R(0)=0,0004983096 $.
Untuk mengana-\breakline lisis jenis kestabilan titik ekuilibrium bebas perceraian dari model (\ref{kompartemen}), substitusi nilai setiap parameter pada Tabel \ref{tbllvariabel} ke dalam persamaan (\ref{uwu}) untuk mendapatkan $\lambda_i<0$ dengan $i = 1, 2, ..., 6$. Karena, $-0,00166<0$, $-0,20166<0$, $0,04264264>0$, $7,700091511>1$, $10,55766496>1$, dan $4,218817709 \times 10^{-9}>0$  maka titik ekuilibrium tersebut stabil asimtotik ketika $R_0<1$. Berdasarkan nilai parameter dan nilai awal itu, grafik untuk setiap subpopulasi terhadap waktu $t$ diperoleh sebagai berikut.
\begin{figure}[h]
  \centering
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/11_M bebas.jpeg}
  \end{minipage}
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/11_V bebas.jpeg}
  \end{minipage}
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/11_Q bebas.jpeg}
  \end{minipage}%
\end{figure}
\begin{figure}[H]
  \centering
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/11_E bebas.jpeg}
  \end{minipage}
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/11_D bebas.jpeg}
  \end{minipage}
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/11_R bebas.jpeg}
  \end{minipage}%
  
  \caption{(a) Grafik subpopulasi \textit{Marriage}, (b) Grafik subpopulasi \textit{Violence}, (c) Grafik subpopulasi \textit{Quarrel}, (d) Grafik subpopulasi \textit{Economy Problem}, (e) Grafik subpopulasi \textit{Divorce}, (f) Grafik subpopulasi \textit{Recovery}}
  \label{solusieo} \vspace{-2.0cm}
\end{figure}
\vspace{+1.2cm}
Selanjutnya, disimulasikan pengaruh variasi $\zeta, \sigma$, dan $\psi$ dengan para-\breakline meter lainnya tetap. Dengan mengambil $\zeta=\sigma=\psi=0,00035$,  diperoleh $R_0= 1,065035325>1$. Karena $R_0>1$ maka perceraian akan meningkat akibat pengaruh individu bercerai. Kemudian, nilai titik ekuilibrium endemik $E_T^*=(1,866559504; 0,001431122832; 0,06720475460; 0,02886420035; 0,01809037588; \\ 0,005801849319)$.

Untuk menganalisis jenis kestabilan titik ekuilibrium endemk dari model (\ref{kompartemen}), agar $\lambda_i<0$ maka substitusi $\zeta=\sigma=\psi=0,00035$ ke persamaan $(3.5) $ sedangkan parameter lainnya tetap sebagaimana pada Tabel \ref{tbllvariabel} ke dalam persamaan (\ref{pkendemk}). Karena, $-1,672967372\times 10^{-21}< 0$ dan $-5,407273372 \times 10^{-32}<0$  maka titik ekuilibrium endemik model (\ref{kompartemen}) $E_T^*$ tidak stabil asimtotik ketika $R_0>1$. Berdasarkan nilai parameter dan nilai awal, grafik untuk setiap subpopulasi terhadap waktu $t$ diperoleh sebagai berikut:
\begin{figure}[H]
  \centering
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/22_Mend.jpeg}
  \end{minipage}
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/22_Vend.jpeg}
  \end{minipage}
   \centering
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/22_Qend.jpeg}
  \end{minipage}%
  \hfill
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/22_Eend.jpeg}
  \end{minipage}
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/22_Dend.jpeg}
  \end{minipage}%
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/22_Rend.jpeg}
  \end{minipage}
  \caption{(a) Grafik subpopulasi \textit{Marriage}, (b) Grafik subpopulasi \textit{Violence}, (c) Grafik subpopulasi \textit{Quarrel}, (d) Grafik subpopulasi \textit{Economy Problem}, (e) Grafik subpopulasi \textit{Divorce}, (f) Grafik subpopulasi \textit{Recovery}} 
  \label{solusie*} \vspace{-2cm}
\end{figure}
Kemudian akan dibandingkan grafik solusi dari model (\ref{kompartemen}) dengan asumsi adanya konseling dan tanpa konseling. Jika diasumsikan tidak ada konseling maka pada model (\ref{kompartemen}) nilai parameter $C=0$. Sebaliknya jika diasumsikan ada konseling maka nilai parameter $C \neq 0$. Dengan menggunakan nilai awal dan nilai parameter yang diasumsikan  $\zeta=\sigma=\psi=0,00035$ dengan parameter lainnya tetap maka pada Gambar 4 diperoleh perbandingan grafik solusi dari model (\ref{kompartemen}) dengan asumsi adanya konseling dan tanpa konseling.
\vspace{-0.5cm}
\begin{figure}[h]
  \centering
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/33_M1 M2.jpeg}
  \end{minipage}
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/33_V1 V2.jpeg}
  \end{minipage}
  \centering
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/33_Q1 Q2.jpeg}
  \end{minipage}%
  \hfill
  \end{figure}
\begin{figure}[H]
  \centering
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/33_E1 E2.jpeg}
  \end{minipage}
    \centering
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/33_D1 D2.jpeg}
  \end{minipage}
  \begin{minipage}{0.3\linewidth}
    \centering
    \includegraphics[scale=0.3]{ARTIKEL JMUA Miya Qarlina Hutagalung/Gambar/33_R1 R2.jpeg}
  \end{minipage}%
  \caption{Perbandingan (a) Grafik subpopulasi \textit{Marriage}, (b) Grafik subpopulasi \textit{Violence}, (c) Grafik subpopulasi \textit{Quarrel}, (d) Grafik subpopulasi \textit{Economy Problem}, (e) Grafik subpopulasi \textit{Divorce}, (f) Grafik subpopulasi \textit{Recovery} dengan asumsi ada konseling dan tanpa konseling} 
  \label{solusiperbandingan}
\end{figure}
\section{Kesimpulan}
Berdasarkan analisis kestabilan sistem dapat disimpulkan bahwa jika $R_0 < 1$, model (\ref{kompartemen}) adalah stabil asimtotik pada titik ekuilibrium bebas perceraian ($E_T^0$). Akibatnya, pengaruh individu bercerai terhadap peningkatan perceraian perlahan-lahan menurun dan akan menghilang dari populasi seiring berjalannya waktu. Namun, jika $R_0 > 1$, model (\ref{kompartemen}) adalah tidak stabil asimtotik pada titik ekuilibrium endemik ada perceraian ($E_T^*$). Akibatnya, pengaruh individu bercerai terhadap pe-\breakline ningkatan perceraian tidak akan stabil dan mungkin tidak bertambah seiring berjalannya waktu.

Hasil simulasi numerik memperlihatkan bahwa adanya pengaruh konseling pada subpopulasi \textit{Divorce}, dimana jumlah individu yang bercerai lebih banyak saat tidak melakukan konseling. Dengan kata lain pada dasarnya konseling memiliki pengaruh terhadap intensitas penurunan jumlah kasus perceraian.

\section{Ucapan Terima kasih}
Terima kasih disampaikan kepada Universitas Andalas yang telah mendanai penelitian ini melalui kontrak Penelitian Skripsi Sarjana (PSS) dengan Nomor kontrak: 258/UN16.19/PT.01.03/PSS/2024.

\begin{thebibliography}{0}
\bibitem{A} Gweryina, R. I., Kaduna, F. S., dan Kura, M. Y, 2021, Qualitative Analysis of A Mathematical Model of Divorce Epidemic with Anti-Divorce Therapy, \emph{Engineering and Applied Science Letters},  \textbf{Volume} : 4(2).

\bibitem{B} Badan Pusat Statistik (2024), "Statistik Indonesia 2024".Tersedia: https://www.bps.go.id/id/publication/2024/02/28/c1bacde03256343b2bf769b0\\ /statistik-indonesia-2024.html. Diakses pada tanggal 23 Mei 2024.

\bibitem{C} Purnamasari, I. A, 2019, Layanan Bimbingan Konseling Keluarga untuk Me-\\ minimalisasi Angka Perceraian, \emph{Irsyad: Jurnal Bimbingan, Penyuluhan, Konseling, dan Psikoterapi Islam}, \textbf{Volume} : 7(1).

\bibitem{D} Gambrah, P. P., Abdul-Rahaman, A. S., dan Adu, A, 2018, Mathematical Model for Minimizing Divorce Through Counseling, \emph{International Journal of Statistics and Applied Mathematics}, \textbf{Volume} : 3(3).

\bibitem{E} Tessema, H., Haruna, I., Osman, S., dan Kassa, E,2022, A Mathematical Model Analysis of Marriage Divorce, \emph{Commun. Math. Biol. Neurosci}, \textbf{Volume} : 2022.

\bibitem{F} O. Diekmann, J. Heesterbeek, dan M. G. Roberts, 2010, The Construction of Next-Generation Matrices for Compartmental Epidemic Models, \emph{Journal of The Royal Society Interface}, \textbf{Volume} : 7.

\bibitem{G} E. Hendricks, O. Jannerup, dan P. H. Sørensen, 2008, Linear Systems Control: Deterministic and Stochastic Methods, \emph{Springer}.

\bibitem{H} Lynch, S, 2007, \textit{Dynamical Systems with Applications using Mathematica}, Boston, Birkhäuser.

\bibitem{I} Bahri S, Pratama EA., 2023, Pemodelan Dinamika Pejudi Muda. \emph{MES: Journal of Mathematics Education and Science}, \textbf{Volume} : 9(1).

\bibitem{J} Diekmann, O., Heesterbeek, J. A. P., dan Metz, J. A. J, 1990, On the Definition and The Computation of The Basic Reproduction Ratio R0 in Models for Infectious Diseases in Heterogeneous Populations, \emph{Journal of Mathematical Biology}, \textbf{Volume} : 28.

\bibitem{K} Bahri, S., Ningsih, S. A., dan Narwen, N, 2023, Kestabilan Model Penularan Penyakit Malaria di Indonesia, \emph{Delta-Pi: Jurnal Matematika dan Pendidikan Matematika}, \textbf{Volume} : 12(2).

\bibitem{L} Diekmann, O., Heesterbeek, J. A. P., dan Roberts, M. G, 2010, The Construction of Next-Generation Matrices for Compartmental Epidemic Models, \emph{Journal of The Royal Society Interface}, \textbf{Volume} : 7(47).
\end{thebibliography}
\end{document}
