{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sakaguchi-Kuramoto model: Adding Phase Shifts to the Kuramoto Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import jax; jax.config.update(\"jax_enable_x64\", True)\n",
    "import jax.numpy as jnp\n",
    "from jax import random\n",
    "from jaxkuramoto import SakaguchiKuramoto, theory, odeint\n",
    "from jaxkuramoto.solver import runge_kutta\n",
    "from jaxkuramoto.distribution import Normal, Uniform\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sakaguchi-Kuramoto model\n",
    "$$\n",
    "\\frac{\\mathrm{d}\\theta_i}{\\mathrm{d}t} = \\omega_i + \\frac{K}{N}\\sum_{j=1}^N \\sin(\\theta_j-\\theta_i+\\alpha)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 0\n",
    "\n",
    "n_oscillator = 10**3\n",
    "alpha = jnp.pi / 3.0\n",
    "loc, scale = 0.0, 1.0\n",
    "dist = Normal(loc, scale)\n",
    "omegas = dist.sample(random.PRNGKey(seed), (n_oscillator,))\n",
    "K = 3.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "t0, t1 = 0.0, 100.0\n",
    "dt = 0.01\n",
    "\n",
    "model = SakaguchiKuramoto(omegas, K, alpha)\n",
    "init_thetas = random.uniform(random.PRNGKey(seed+1), shape=(n_oscillator,), maxval=2*jnp.pi)\n",
    "sol = odeint(model.vector_fn, runge_kutta, t0, t1, dt, init_thetas, observable_fn=model.orderparameter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x107cb9ca0>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0jklEQVR4nO3dd3xV9fnA8c+TCYSQEAgJEDYJYWMIW/bGgbVaxb2q1lFXXbXaWm1rbbXW/UPFvRUFKQiIIHvPhDBCSCBkh+yQeb+/P+5NyL6XwM183q9XXnLOPefc5x7vPc93ne8RYwxKKaVUXVwaOwCllFJNnyYLpZRSdmmyUEopZZcmC6WUUnZpslBKKWWXJgullFJ2OS1ZiMgiEUkRkYhaXhcReVVEokVkv4iEOSsWpZRS58eZNYsPgDl1vD4XCLb93Qm85cRYlFJKnQenJQtjzHrgdB2bzAc+MlZbAV8R6eqseJRSStWfWyO+d3fgZIXleNu6xKobisidWGsfeHl5jQwNDW2QAJVSqqXYtWtXmjHGv777N2aykBrW1Tj3iDFmIbAQIDw83OzcudOZcSmlVIsjInHns39jjoaKB3pUWA4CEhopFqWUUnVozGSxFLjJNipqLJBljKnWBKWUUqrxOa0ZSkQ+B6YAnUUkHvgz4A5gjHkbWA7MA6KBfOBWZ8WilFLq/DgtWRhjFth53QD3Ouv9lVJKXTh6B7dSSim7NFkopZSyS5OFUkopuzRZKKWUskuThVJKKbs0WSillLJLk4VSSim7NFkopZSyS5OFUkopuzRZKKWUskuThVJKKbs0WSillLJLk4VSSim7NFkopZSyS5OFUkopuzRZKKWUskuThVJKKbs0WSillLJLk4VSSim7NFkopZSyS5OFUkopuzRZKKWUskuThVJKKbs0WSillLJLk4VSSim7NFkopZSyS5OFUkopuzRZKKWUskuThVJKKbs0WSillLJLk4VSSim7NFkopZSyS5OFUkopuzRZKKWUskuThVJKKbucmixEZI6IHBaRaBF5oobXfUTkBxHZJyKRInKrM+NRSilVP05LFiLiCrwBzAUGAQtEZFCVze4FDhpjhgNTgJdExMNZMSmllKofZ9YsRgPRxpgYY0wR8AUwv8o2BvAWEQHaA6eBEifGpJRSqh6cmSy6AycrLMfb1lX0OjAQSAAOAA8YYyxVDyQid4rIThHZmZqa6qx4lVJK1cKZyUJqWGeqLM8G9gLdgBHA6yLSodpOxiw0xoQbY8L9/f0vdJxKKaXscGayiAd6VFgOwlqDqOhWYLGxigaOA6FOjEkppVQ9ODNZ7ACCRaSPrdP6WmBplW1OANMBRCQAGADEODEmpZRS9eDmrAMbY0pE5D5gJeAKLDLGRIrI3bbX3waeAz4QkQNYm60eN8akOSsmpZRS9eO0ZAFgjFkOLK+y7u0K/04AZjkzBqWUUudP7+BWSilllyYLpZRSdmmyUEopZZcmC6WUUnZpslBKKWWXJgullFJ2abJQSilllyYLpZRSdmmyUEopZZcmC6WUUnZpslBKKWWXJgullFJ2abJQSilllyYLpZRSdmmyUEopZZcmC6WUUnZpslBKKWWXJgullFJ2abJQSillV53JQkRcRGR8QwWjlFKqaaozWRhjLMBLDRSLUkqpJsqRZqhVIvJrERGnR6OUUqpJcnNgm4cBL6BURM4AAhhjTAenRqaUUqrJsJssjDHeDRGIUkqppstuM5RY3SAiT9uWe4jIaOeHppRSqqlwpM/iTWAccJ1tORd4w2kRKaWUanIc6bMYY4wJE5E9AMaYDBHxcHJcSimlmhBHahbFIuIKGAAR8QcsTo1KKaVUk+JIsngV+A7oIiJ/AzYC/3BqVEoppZoUR0ZDfSoiu4DpWIfNXmGMiXJ6ZEoppZoMu8lCRD42xtwIHKphnVJKqVbAkWaowRUXbP0XI50TjlJKqaao1mQhIk+KSA4wTESyRSTHtpwCLGmwCJVyQEFxKacyzzR2GEq1WLUmC2PMP2x3b//LGNPBGONt++tkjHmyAWNUyq57Pt3NhBd+ZvOxtMYORakWyZFmqKfqewe3iMwRkcMiEi0iT9SyzRQR2SsikSLyyznErhQAqTmF/HwoBYCvd8Y3cjSqOckuKMYY09hhNAuOJIs3qMcd3La+jTeAucAgYIGIDKqyjS/WO8QvN8YMBq52OHKlbHafyACgu29bthxLb7D3zcov5q6Pd/L1zpMN9p7qwvl650mG/WUV172zjcW74ykubZ63j+UVlvDqmqMcTc5x6vs4kizGGGPuBQrAegc34Mgd3KOBaGNMjDGmCPgCmF9lm+uAxcaYE7ZjpzgcuWqxfoxI4p5Pd7HHlgTs2R2XgYerC9eN6UlSdgGn84ouWCxpuYWk5BTU+Non2+JYGZnMo9/sJ+tM8QV7T9UwFu8+BUBkQhYPf7WPV3460mDvXVhSWutrGXlFPLfsIOuPpDp0rA82x/Ly6iM89NXeCxRdzZx5B3d3oGKRK962rqIQoKOIrBORXSJyU00HEpE7RWSniOxMTXXsBDZlqTmFrD2c0uyrvxaL4UB8FgXFtX/xY9PyzulCGp2Sw++/2MPyA0nc+N52olNy7e6zKy6DoUE+DO3uA8ChpGyH368upRbDb97ewqQX15KaU1jttc+2naCdhysAPx9KviDvqawOxGcx6m8/ceWbm8i4gMm/TE5BMTtiT3P35H7sfnom84YGsmhjbLX/zxdaqcXw8Fd7CX36R97fdLzGbZ5YvJ/3Nh7njg93su9kpt1jboq29tNFnMomJtX+76W+6nsH998d2K+mhyVVvTq6YR2GewkwG3haREKq7WTMQmNMuDEm3N/f34G3btpuXrSdW9/fwbL9ifU+Rl5hCTe8u435r2906ILqDE8vieCy1zfyqzc3k19UUr7+RHo+x9PyeOybfUz59zqmv7SOzHz7P/jTeUU89OU+vDxcWXrfBDzdXLjz453kFZbUuk9BcSn747MY2asjoV2ts+kfTqpfdbyk1EJMai4ltuaIn6KSiUnLo6DYwvd7TpVvdygpm/+sPsKpzDO8dPVwAjp4sjLCuckir7CE/+1PrHSez5Uxhvs+280jX+3DYnG8oHKmqJTX1hzl5kXb+ffKw+dUyDmWmstT3x3gpVWHeWnVYd7dEENRif2y5j9WRJFTUMyBU1k897+DDr+fozZFp1FiMUwL7YKbqwuPzg6lqNTCq2uOUlxqOe+C3Fc7TzLir6v4v1+OVVq/aONxFu8+RXsPN15adaRaQSsqMZuVkcncPK4X/t6e3PL+do7U0bxkjCEqMZvJIdbr4oqIpPOKuy52k4Ux5lPgMaxTfCRivYP7aweOHQ/0qLAcBCTUsM2Pxpg8Y0wasB4Y7kjgzVVydgEHE60l3y92nHBon30nM3l51eHyUo/FYvjD1/vYGJ1GZEI2v/m/LQ5XWWuyKy6DG9/bxvbjpx3eJzm7gM+3n2BED18OJWXzxtpoAJKyCpj36gam/nsdX+2M55KhXUnLLeL9TbF1Hi86JYcZL//C4aQc/vnrYQwL8uW1BRcRk5rHp9viat1vf3wWRaUWRvX2w7+9J35eHnX+uGpTUFzKZa9vYtpLv7Dgna0UlVh4d0MMQR3b0s/fi19s5/ep7w4w55UNvL42muE9fJk1OJCZgwL45UhqnTWs8/XsD5Hc+9lunvj2AGAtGW+LSaf0HC76e05msmx/It/ujufHyLMXlY+3xHLHhztIz61eqk7KKmDBO1t5+acjHEvN5fW10Xy7+1S17aoyxvDamqPMeWU9X+44yRtro3l9bTTP/y+K/9hp7snML2JLTDp3TerH9WN6sXRvgsPDoi0WQ9aZYlKyC9hbR6l87aFUvNu4EdbTF4A+nb1YMLoHH2+NY+DTP3LlW5vLCw3nKioxm6e+O0BmfjH/XXOUnAJrzTouPY+XVh9mxsAuvHXDSHILS9hwtPLovY+3xtHW3ZWHZobw2W/H4O7qwg3vbqu1sJWUXUBGfjHTB3ZhWJBP+UAPZ3CkZgGQDGwANgNtRSTMgX12AMEi0sc2S+21wNIq2ywBJoqIm4i0A8YALXoqkbIO2KkD/NlyLN1utTcuPY/r3tnKqz9Hc+3CLfx9eRSXv7GRFRFJPDVvIKsfnkw7D1duWrSdL7Zbk8+6wyn89YeDDl28jDH8eWkEG46m8Yev9zl88dl8LA2Lgb//aiiXDO3K+5tiScoq4G/Lo8gtLGHGwADun9af16+7iOmhXfh4a1yd8fzLVmJdct8EZg0OBGB8/86E9+pY3rZckx2x1gQX3qsjIkJIQHsO1aNm8dGWWKISs/lNeBA7YjO4ZuEWdsRmcPvFfZgc0oXtsaf56WAyn247wfVjerLigYl8c/c4XF2E2YMDOVNcyrrD9hN2bmEJ++MzeWZJBCsjHSsFlpRa+NFWYly6L4GIU1nctGg71yzcyoNf7uXHiCSHLqafbI2jjbsLXbw9ee3naCwWw+4TGTy9JJKfolJ4Z0PlZpGUbGviP5yUw1vXj2TDY1MZ2LUD76yPsVvyXrz7FC+tPsKswYFs/eN0Ip+dw9Hn53LJ0K58vv1End+zLcfSMQYmhXTmt5P6AjD5xbU8sySi1vdNzDrDI1/tY+Tzqxn+7CpG/30NV7yxqcZClDGGtYdTmBTij5vr2Uvg43NCuWtyXyaF+LPnRCbrj9avAPbZthO4ubjw3s3h5BeVlv+/e+Wno7iI8NwVQwjv3REPV5fy7y9YE92qyGSmhXbBt50HvTp5seiWUaTlFvLqmuga3+tggrXgObBrB6YM6MKeExlOabYDxx5+9BywH2tz1Eu2v3/b288YUwLcB6zEmgC+MsZEisjdInK3bZso4Efb8bcD7xpjIur5WZqFzcfS8GnrzmNzQrEYWBFRc1NUYUkpR5NzeOyb/biI8PwVQ4hJy2Ph+hgAHpkZwh0T+9Cnsxc/PTyZkb068spP1lLM/Z/tYdGm43y6rfaai8ViyCko5pcjqUScymbKAH9OnM5n3eGaSyZVmw62HEvHt507oYHePDY7FIsxzPzPL/ywL4EHpgfz7s3hPDJrACLCHRP7cjqvqNaLfmFJKRuOpnHJsK4M7Fr5ab2zBwdyKCmH+Iz8GvfdEXua4C7t6ehlHXMxIMCbI0k51S4qdTV9ZOQV8cbaY0wO8efFq4Zz3Zie7DmRyeBuHVgwuidTBvhTVGLhjo920qtTO56+dBADu3bA3XahGdu3E1192vDgl3t4cvH+8lE1Px9KZsILP7Nkr/VzP7fsIEP+vJLLX9/ER1viuP+zPcSm5dUaV5ldcRlkF5Twz18PxdvTjUtf28ieE5mEBnrzw74E7v5kF9e9s5X8ohK2xaQz5V9ruefTXZUuyBl5RSzbn8ivw4J4cl4oUYnZPPTVXm56bzvdfNoQ3qsjKyISK5239zYdJ/tMMYvvGc+cIYGICDeN68Xh5Bz2x2dVizMrv5i/LI1kyd5T/G15FCN7deS1ay+ic3tP2nq44ubqwqzBAWTmFxNxqvr+ZTZGp9He041hQb50923Ls/MHM7BrBz7aEseSvVUbJ6zfnzs+3Mn/DiQwdUAXnpwbyoMzgnF3FT7aEltt+8iEbFJyCpk6oEul9d5t3Hly7kDeuiGMNu4u/FIl+afmFJJdUHf/m8ViWHUwickh/kwL7UKvTu34fu8p4jPyWbovgQWje9LVpy1t3F0ZFuRTqTa/52QmabmFzBocUL5uSHcfrhnVg4+2xLI1pvpIv7KC0YBAb6YO8MdiqHeSs8eRmsVvgH7GmCnGmKm2v2mOHNwYs9wYE2KM6WeM+Ztt3dvGmLcrbPMvY8wgY8wQY8wr9foUzciWmHTG9vVjYNcODAjwZmmFL39JqYUVBxLZfCyN2f9Zz8z/rGfb8dM8MiuEG8b24uu7xvHt78ax7P6J3D89GBFrt1Abd1fundqPpOwCfv/5HnJsbfzvbYipcTjgtph0Zrz8C0P/sopb3t9BN582vHl9GJ3be/BVDcNAY1JzGfHXVfzxuwOVPseYPn64uAg9O7Vj4Y3hjO7tx+NzQnlgenCl/cf29WN4kA+vrjlaY//DtpjT5BeVMi20S7XXpoZa22J/qaGEWGox7IrNYFQfv/J1IYHe5BVVvpv78+0nCH16BU9ViL+i/645Sm5hCU/OCwXg+flDWHb/xSy+Zzxt3F2Z0L8zPf3aAfDc/CG0cXettL+7qwsf3jaa2YMD+Xz7SZ5bdpCSUgt//eEgpzLP8MySSL7bE897G49zydCu/PfaESz//UQ83Fx49odIu6X0NYdScHcVLhnWjdsn9gFgUoh/ee3mT5cMJC49nwkv/Mz1724j7nQ+yw8klTdzGmN4ceVhikos3Dy+N5cN68akEH+W7E1gRA9fvrxrHFeGBRGXnl9+8bFYDEv2JDA5xL9SAr9kWFc83Vz4dnf1+1neXBfNB5tjeeCLvWTmF/Hc/CG4uFTuuhzfrzMAG+q4oG0+Zv1ulSXj68f0Ysm9Exja3YcXfzxUqYZ6pqiU3360i8iEbF5bEMbL14zgrsn9eHBGCNeP6cWGo2nV+nnKCkRl7fxVebq5MqZPJzZGn20i2hydxvgX1jD7P+vLm5Vqsjc+k+TsQmYPCUBEuGJEdzYfS+e5ZQcR4PaL+5RvO7J3RyITzg4QWRWZhJuLMKVKEntk1gC6+rbh+ne3VaqJgLV/rrtvWzq0cWd4kC+dvDyc1hTlSLKIAHyd8u6tzMnT+Zw8fab8B3P5iG7sjMvgRLq11Pzqz9H87tPdXPfONrILSvjHlUP58LbR3Dy+NwDhvf0Y2cuvxmNPCelCn85erD2cSq9O7Xj3pnASsgr4bk/l0vwP+xK48b3tGODR2QP47cQ+LLwpnHYeblwZFsSaqBTSqrRdf783gfyiUj7bdoJTmWc4lXmGk6fPMLZvp/JtJoX4894to/jdlH7VLhAiwjOXDSIpu4DXfq5enf75UAqebi6M69u52mv9/NvT3bdtjU08BxOyySksYVTvjuXrBgRYO7nL+i2yC4p5ftlBLAY+3XaiWkk+I6+IL3ac4KqwIEIDrRdFFxdhSHcfPN2sScHVRVh8z3hWPjiJSbVcYEICvPnvtRdx56S+fLQljqve3kJsej4Pzwwh60wxD325j9BAb16+ZjjzR3RnULcOPDgjmLWHU7n/8z21lliNMaw+mMzYvp1o7+nG/dOC+eT2Mbx1fRgiQnhvP+6Y2JcFo3uQkV/MhP6d2fvMLMb08eOlVUfIOlPMOxti+Hz7CX47sQ8hAd64ubrw4a2j2PvMTD65Yww9/Noxc1AAIpQ3je2IPU1SdgHzL6o8gLFDG3dmDQ5k6b6EasM/1x5OYVTvjjw5N5SPbhvDoG6Va4kA/t6eDO7WgfVHrBfi4lJLpc72pKwCjqflMa5fp0r7ubgIf5w3kISsAhZtOk5xqYXTeUXc+sF2NhxN5cVfD2PmoIBK+0wf2IXCEku1e2/WHk5leJAP/t6eNZ5zgIv7d+ZYah5JWQUUlpTy6Df7KbUYErMKKg12qGql7YI/LdQay5Vh3XERYWVkMr+6qDvdfNuWbxvWsyPFpYaIU1kYY1gZmcS4fp3waete6Zid23uy7P6J+Lf35NU1Ryu9djgphwGB3uXnaHz/zmyNSXfKSEtHksU/gD0islJElpb9XfBIWoGyL23ZD+FXth/i93tPkZlfxKKNxxkW5MOzlw/m+3smsGB0TyaH+JfXIOri4iI8OnsAPfza8qdLBjEttAsjevjy/LKDJGZZS9lHknN4+Ku9DO/hw/f3TODeqf156pJBDLENOb16ZBAlFsN172zlL0sjy5tufoxIpLtvW0Tgqx0n2Wr7HBWThT0je/lxZVh33t90vFI/TVn78fh+nWjr4VptPxFhygB/NkenVWtK+vlQCiIwMfjsBTzE9sMpKyF/vTOevKJSFt44EhGqJc/Pd5ygoNjCbRVKfDXp3N6z/EdZl8fnhDJ3SCB7T2Yypo8f90/rz/VjeuLT1p2XfjO8PAEB3DahD4/OHsCKiCQufXVjjcMk391wnONpeVwxwvpdcXURLg7ujJdn5Qmj/3zZYD7/7VgW3TIKn7buPH3pIDLyi7jm/7bwwopDzBsayB/nDSzfXkTwbXf2dil/b09G9fJj6b4Eaw03IgkPNxem11DbuzKsO5n5xaytUIJNySngSHIu0wcGcNfkflwcXD3xl5kywJ/tsad54tv9jHxuNTe8t608YWw7Xvt3a1y/TswYGMCLPx4m7K+rCXtuNTtiM3jlmhH8ZlSPatuP7uNHW3fXSgWNjLwi9pzIqFZ6r2p8f+v7b4pO4/s9pziVeYYPbh3NoK4d+HpXPFlnivli+4lKnejGWPscKl7we3Xy4p2bRnLD2J48MTe00nuE9bQWcnafyCA6JZfY9Pzy/rqqfNq6c+3oHmyMTuPkaWvhMjO/iCMpOQwL8jn7mXt3JDm7kPiMyn1Yzyw5/9Z9R5LFh8A/gRc422fx0nm/cyv0U1QyAR08Ce7SHoBuvm0Z1bsjy/Yn8M6GGHILS3jxqmHcPL43PTu1O+fjzxvalQ2PTWPmoABcXIRXrhlBQYmFl1YdodRieOyb/Xi3ceftG0bi08692v7BAd5cNrwbx1Lz+GBzLJ9sjSM6JZcjybncMbEPk4L9+XLHSdYcSsbPy6O8FO+o+6b2p6jUwqIK48tj0vKIS89nag0XpTKTQ/zJKyqt1mb786FkRvTwpXP7syXEDm3cCerYlt1xGZSUWvhg83FG9urIrMGBjOvbiaX7EspLXcWlFj7aHMfF/Ts7lAgc4eoivHFdGEvuncCHt41GbP1Nu/40g8HdfCpt6+Ii3Du1P1/eOZaSUgtXvb258gU4u4B/rTzMrEEBXBlW9Ralytq4uzKuXydcbbW6Id19+HVYEIeScgju4s2LVw23W+i4eXxvYlLzeODLvXy35xSTQ/yrJSWAif074+/tWakPqqwgNL6f/QLE1SN70M7DlW92xdPRy4PNx9LLazRbjqXj3catWt9VmT9fNohOXh50823L76cH89VdY5k/ouZz4+lmbUKseE/Tsv0JWAzMGBhQ4z5lBgZ2wM/Lgw1HU/m/X2IY0r0DE4M7c3V4EPvjswh/fjVPLD7AlW9uYlec9ebRoym5HE/LY3aVC/600ACev2IondpXrsn4e3vS068du+Iyyj//zDriujq8B8LZUZQ/RaVgjLUWVCa8t7XloWJz1daY9Dr7Lx3lSLJIM8a8aoxZa4z5pezvvN+5lckpKGbdkVTmDula6Ud76bBuHEnO5Y21x7hseLfyppALoXdnL24c24vFu+N58Mu97D2Zaf2xta+9+v3agos49vd5jO3rx5vrjvHWumO4CFwytCsLRlvvkF5+IIlLh3Wt1txkT1//9swb0pVPtsSVDwVcbrvXpK4f76QQf3zaulfqT0nKKmBffFaN+80ZHMhPUSmM+fsaTp4+w522ETWXD+/G8bQ8Ik5ZR5As259AUnYBt07ofU6fwx4XF2F4D9/yvg0RqTTqpqrw3n4sf2AiAwK9+e1HO/njdwfIzC/i461xFFssPHXJQIdql1U9N38I/712BJ/fOZb2NVz0q5o3NJAHpgfzv/2J5BQUc/fkvjVu5+bqwoyBAWyJSa9QIziNt6dbtYRYk96dvfjl0als++N0fn5kCl28PVm85xQlpRZ+ikpmUrB/edKrqodfO7Y/NYOVD03i4ZkhtTbLlpkywJ/4jDNEp+Ty+8/38PSSSIb38GVI97p/Zy4u1hrt93sTiEnL4+7J/RAR5o/ojpeHK8bAK9eMwLedB+9ttA46WbY/ERGYNajuRFRRWE9fdp/IZGWkteAT6NOm1m27+7ZlWmgX3t8Uy43vbePxb/cTGujNyF6Vm2G927ixI9aawD7YdJzr3tlKUMe2tR3WYfa/QbBLRP6BddhrefuBMWb3eb97K/JTVDJFJRYuG9610vp5Q7vy75WHcXERHps94IK/7z1T+vH9nlP8sC+B+SO6cfnwbg7td/+0YK5/dxvf7o7nxrG96NKhDdMHdqGfvxcZ+cXcPblfveK5d2p/VkYmccv7O7huTE8Wro9hQv9Oldpyq2rj7sqvw4L4eKv1Dlt/b09WR1lvgps9uPoP875p/cktLKGoxEJ4b7/yH+/cIV15ekkES/edIqewmL8sPUhooHe1UTGNwbedBx/eOpqXVx/hyx0n2XMikxPpecwYGECvTl71OmZbD9daS901EREemhnCvKFdsRhTa+ke4KKevny+/QQxaXl0923Lz1EpjOnbqdaLfFUV+wsuH96ND7fEsiIiibTcomq/kaocfQ+wJguA33+xl6jEbK4J78HDs0IcSr73TwtmU3QaAwI7MMdWW/Dz8uDHBydRWGKhf5f27Iw7zTe74skrLOGHfQmM69uJLh1qv+BXNbJXR77fm0BqTiGPzbH/+39wRghbjm0hKjGHG8f24vaL+1T6LC4uQnivjmyMTuWDTcf5yw8HmTUogH//Zjg+jzkcVo3EXkeIiKytYbVxdETUhRYeHm527tzZGG99Xm7/YAdRidlsfHxatRL56bwiLMZUak65kNJzCzmSnMvoPn7n9ENbcSARA8y1DZsE681rIlRqez9XqyKTuO+zPRSVWggN9Oadm8Lp4Vd3s5v1pr31PDZnAPdM6c9tH+wgJjWXtX+Yck6l7js+3Mm6wykYrDdiLbp5VL2a/JxpxYFEfveptSz27e/G2S09N4ajyTnM/M96XrxqGP/bn8j6o6l8dNvoSv1Hjtp3MpP5b2wCoHN7DzY+Pq3aiLPzMfPlXziaksvwHr5897vx51QjNsbU+f3aFpPONQu3ctuEPizadJwXrhzKtaN7Onz8Y6m5TH/J2lCz+YlpdRaayuQWluDuKrX+Bv+3P5F7P7N+f2YOCuDN68Nwd3VBRHYZY8IdDq4KR57BPbW+B1dWmflFrD+ayi3je9f4RfXzcmRexvrr1N6TcfVIRHOHVi/hXYgf8azBgax+eBIxqXlMDO5cZxNNmf5dvBnTx4/Pt5/gtgl92HIsnavDg865eeZPlwykoLiUrj5tePqyQXRoU73vprHNHdqV924Op7jU0iQTBVhHqXVs587fl0eRmV/M05cOqleiABgW5EMPv7acPH2Gm8b1vqCJAuD+6cF8vCWW568Yes5Np/a+X6N6+9Hdty2LNh3H082FOUNq7qCuTT//9rxzUzjtPd0cShSA3SbFeUMD+cOsEHILS3loZnD5EOTz5UgzFCJyCdbHq5bXr4wxf70gEbQCP+xPpLjUcMVFjjcJtHS9Onmdc/PKgtE9efDLvfx3zVHOFJdW6thzVO/OXnxyx5hz3q+hTbfTAdvYXFyEicH+LN2XQO9O7bhpXK96H0tEeG1BGJuPpXHbhLpHpdXH5cMdb349Vy4uwu+m9ONP30dwy/jelUaYOarqkN/zJSLcNy3Y/obnyG6yEJG3gXbAVOBd4Cqsd1srB+QXlfDh5lhCA70ZVEcbsLJvzpBAOi3z4K11x/Bt517rPQ+qYdw1uS/xGfk8MXfgeZdeR/TwZUQP3wsTWAO7YWwvZg0OwN9JzchNhSP/h8cbY24CMowxz2J9EFL1Qc2qmuJSC9e/u42Y1FwenT2gXiNa1Flt3F15+8aRTAzuzAtXDrvgzRXq3Azu5sPieyYwuk/TbCprSF2827T437cjzVBlT37JF5FuQDpw4euKLdB3u0+x50QmL/9meJNvVmguRvX24+Pbm34zklItjSPJ4gfb40//BezG+kyKd5wZVEuxeE88/bu0L79TWymlmqs6k4WIuABrjDGZwLcisgxoY4ypfcpIBVhHQG0/fpp7p/Zv8dVTpVTLV2efhTHGQoWpPYwxhZooHBOZkI3FnNv8SUop1VQ50sG9SkR+LVo8PidRtqfhhV6gOYeUUqoxOdJn8TDgBZSISAHWZ2sbY4yOA63DoaQc/L0965yHSSmlmgtH7uDWonE9HErK1lqFUqrFcPQO7o5AMJXv4F7vrKCau5JSC0eTc8sfWqSUUs2dI3dw3wE8AAQBe4GxwBagUSYSbA5i0/MpLLGc8/MelFKqqXKkg/sBYBQQZ5tU8CLAOU8EbyEOJdk6t7tqslBKtQyOJIsCY0wBgIh4GmMOARf+wQstyKHEHFxdhP62J+IppVRz50ifRbztDu7vgdUikgEkODOo5u5QUjb9/L3O65kPSinVlDgyGupXtn/+xfYgJB/gR6dG1cxFJeYQVuFRh0op1dw5NK+wiISJyO+BYUC8MabIuWE1X1n5xZzKPMNA7a9QSrUgdpOFiDwDfAh0AjoD74vIn5wdWHMVZevcruv5xUop1dw40mexALioQif3C1hnn33emYE1V2XTfOiDjpRSLYkjzVCxVLgZD/AEjjklmhbgUGIOfl4edPHWaT6UUi2HIzWLQiBSRFZjfZbFTGCjiLwKYIz5vRPja3aikrIZ2NVbpyVXSrUojiSL72x/ZdY5J5Tmr6TUwuGkHG4cW/+H1yulVFPkyNDZDxsikJYgNj2PwhKLdm4rpVoch4bOKsccTMwBdCSUUqrl0WRxAUUlZuPmIvTr4tXYoSil1AVVZ7IQEVcR+VdDBdPcHYjPIiTAW6f5UEq1OPaewV0KjNRHqtpnsRj2x2cyoqdvY4eilFIXnCOjofYAS0TkayCvbKUxZrHTomqGYtPzyC4oYXiQT2OHopRSF5wjfRZ+QDrWhx1dZvu71JGDi8gcETksItEi8kQd240SkVIRucqR4zZFu+IyABjew7dxA1FKKSdwZOjsrfU5sIi4Am9gvYkvHtghIkuNMQdr2O6fwMr6vE9T8cuRVPy9PfXpeEqpFsmRiQRDRGSNiETYloc5OJHgaCDaGBNjm6X2C2B+DdvdD3wLpJxD3E1KSk4Bqw4mM3twgN65rZRqkRxphnoHeBIoBjDG7AeudWC/7sDJCsvxtnXlRKQ78Cvg7boOJCJ3ishOEdmZmtr0nui6aGMsJaUW7ri4b2OHopRSTuFIsmhnjNleZV2JA/vVVMQ2VZZfAR63jbqqlTFmoTEm3BgT7u/v78BbN5yiEgufbotj7pCu9O6s91copVomR0ZDpYlIP2wXelsndKID+8UDPSosB1H9cazhwBe2ppvOwDwRKTHGfO/A8ZuEiIQscgpKuHRY18YORSmlnMaRZHEvsBAIFZFTwHHgBgf22wEEi0gf4BTWpqvrKm5gjOlT9m8R+QBY1pwSBcDO2NMAjOytj1FVSrVcjoyGigFmiIgX4GKMyXHkwMaYEhG5D+soJ1dgkTEmUkTutr1eZz9Fc7EzNoNendrRxbuN/Y2VUqqZqjVZiMjDtawHwBjzsr2DG2OWA8urrKsxSRhjbrF3vKbGGMOuuAwmD2ha/ShKKXWh1VWzKLthYAAwClhqW74MWO/MoJqL2PR80vOKCO/l19ihKKWUU9WaLIwxzwKIyCogrKz5SUT+AnzdINE1cWX9FeHaX6GUauEcGTrbEyiqsFwE9HZKNM3MrrgMOrRxo79/+8YORSmlnMqR0VAfA9tF5Dusw2d/BejT84AdsacJ7+2Hi4veta2UatnqTBa2qck/AlYAE22rbzXG7HF2YE1dem4hx1LzuGpkD/sbK6VUM1dnsjDGGBH53hgzEtjdQDE1Cztts8yO0v4KpVQr4EifxVYRGeX0SJqZHcdP4+HmwlB9foVSqhVwpM9iKnC3iMRiffiRYK10DHNmYE3d9tjTjAjy1UeoKqVaBUeSxVynR9HMpOUWsj8+i4dnhjR2KEop1SDsNkMZY+IAX84+Jc/Xtq7VWn/EOk361AFdGjkSpZRqGI48/OgB4FOgi+3vExG539mBNWVrD6fSub0ng7t1aOxQlFKqQTjSDHU7MMYYkwcgIv8EtgCvOTOwpupMUSnrDqUwZ0ig3l+hlGo1HBkNJUDFhxOVUvODjVqFHyMTySks4cqwoMYORSmlGowjNYv3gW22O7gBrgDec1pETdzXO+Pp6deOMX108kClVOvhyPMsXhaRdcDFWGsUrfYO7sSsM2w+ls7DM0O0CUop1ao4UrPAGLMbvYOb1QeTAbhEH6GqlGplHOmzUDarIpPp6+9FP51lVinVymiycFDWmWK2xqQza1BgY4eilFINTpOFg9YdTqHEYpg1OKCxQ1FKqQanycJBKyOT8Pf2ZESQb2OHopRSDU6ThQPyCkv4+VAKcwbrjXhKqdZJk4UDfopKpqDYwmXDuzV2KEop1Sg0WThgTVQKndt7EN5LH3SklGqdNFnYYbEYNkanMTHYX5uglFKtliYLOyISsjidV8SkkM6NHYpSSjUaTRZ2lD27YmKwfyNHopRSjUeThR3rj6QxpHsHOrf3bOxQlFKq0WiyqENOQTG7T2QwSWsVSqlWTpNFHTYfS6fEYpgUoslCKdW6abKowy9HUvHycCWspw6ZVUq1bposamGMYf2RVMb164yHm54mpVTrplfBWhxLzSU+4wyTB2gTlFJKabKoxYoDSQDMHKizzCqllCaLWqyISCKspy+BPm0aOxSllGp0Tk0WIjJHRA6LSLSIPFHD69eLyH7b32YRGe7MeBx1Ij2fg4nZzB2ij09VSilwYrIQEVfgDWAuMAhYICKDqmx2HJhsjBkGPAcsdFY852JFRCIAc4boU/GUUgqcW7MYDUQbY2KMMUXAF8D8ihsYYzYbYzJsi1uBICfG47Af9icwLMiHHn7tGjsUpZRqEpyZLLoDJyssx9vW1eZ2YEVNL4jInSKyU0R2pqamXsAQqzuanEPEqWyuGFFXqEop1bo4M1nUNJ+3qXFDkalYk8XjNb1ujFlojAk3xoT7+zt3KOviPadwdRF90JFSSlXg5sRjxwM9KiwHAQlVNxKRYcC7wFxjTLoT47Er60wxn2yNY8bALvh768SBSilVxpk1ix1AsIj0EREP4FpgacUNRKQnsBi40RhzxImxOOT9TcfJKSjh/mnBjR2KUko1KU6rWRhjSkTkPmAl4AosMsZEisjdttffBp4BOgFvighAiTEm3Fkx1SXrTDHvbTzOrEEBDOnu0xghKKVUk+XMZiiMMcuB5VXWvV3h33cAdzgzBke9uS6anIISHpihtQqllKpK7+AGolNyeH9jLFeGdWdwN61VKKVUVa0+WaTlFnLnR7vwbuPG43NCGzscpZRqkpzaDNXUFZdauO2DHSRkneHj28cQ0EHngVJKqZq06mTxwaZY9sdn8eb1YYzq7dfY4SilVJPVapuhikosvLsxhgn9OzFvqE4YqJRSdWm1yWJlZBLJ2YXcMbFvY4eilFJNXqtNFv/bn0hAB08mB+uT8JRSyp5WmSzOFJWy7kgKswcH4uJS0xRWSimlKmqVyeKXI6kUFFuYM1ifV6GUUo5olcliZWQSvu3cGd1HR0AppZQjWl2yKC61sCYqmemhAbi5trqPr5RS9dLqrpY7YzPILihh5qAujR2KUko1G60uWayJSsbD1YWJOgpKKaUc1qqShTGG1VHJjO/fCS/PVn3zulJKnZNWlSz2xWcRl56vo6CUUuoctapk8c2uk7Rxd2HeMJ3eQymlzkWrSRZpuYV8u+sUlwztRoc27o0djlJKNSutIlmUlFp4btlBCktKuWdqv8YORymlmp0W3cubmV/ECysOsTIyiYz8Yh6ZGUI///aNHZZSSjU7LTZZFJaUcsv7O4hMyOKyYd2YMySQmYMCGjsspZRqllpssnj2h4PsPZnJ2zeEMWeIdmgrpdT5aJF9FuuPpPLZthPcNamvJgqllLoAWlyyyCss4Y/fHaCvvxcPzQxp7HCUUqpFaHHNUC+vPkJ8xhm+umscbdxdGzscpZRqEVpUzWLfyUze33Sc68b01OnHlVLqAmoxyeJ0XhEPfLEHf29Pnpgb2tjhKKVUi9IimqF+jEjimSURZOYX8/mdY/QObaWUusCadbJIySngz0siWRGRxMCuHXjv5lEMDfJp7LCUUqrFabbJYt3hFB75ah85hSU8OnsAd07qi7s++U4ppZyiWSaLL3ec4MnFBwgJ8OaLO8cSHODd2CEppVSL1uySxem8Ih7/9gCTQvz5vxtG0tZDh8cqpZSzNbt2m1OZZ5gywJ93bwrXRKGUUg2k2SWLtu6uvLbgIjzcml3oSinVbDW7K24///Z469BYpZRqUM0uWYg0dgRKKdX6ODVZiMgcETksItEi8kQNr4uIvGp7fb+IhDkzHqWUUvXjtGQhIq7AG8BcYBCwQEQGVdlsLhBs+7sTeMtZ8SillKo/Z9YsRgPRxpgYY0wR8AUwv8o284GPjNVWwFdE9AEUSinVxDjzPovuwMkKy/HAGAe26Q4kVtxIRO7EWvMAKBSRiAsbarPVGUhr7CCaCD0XZ+m5OEvPxVkDzmdnZyaLmrqiTT22wRizEFgIICI7jTHh5x9e86fn4iw9F2fpuThLz8VZIrLzfPZ3ZjNUPNCjwnIQkFCPbZRSSjUyZyaLHUCwiPQREQ/gWmBplW2WAjfZRkWNBbKMMYlVD6SUUqpxOa0ZyhhTIiL3ASsBV2CRMSZSRO62vf42sByYB0QD+cCtDhx6oZNCbo70XJyl5+IsPRdn6bk467zOhRhTrYtAKaWUqqTZ3cGtlFKq4WmyUEopZVezShb2pg9pyUSkh4isFZEoEYkUkQds6/1EZLWIHLX9t2Njx9oQRMRVRPaIyDLbcms9D74i8o2IHLJ9N8a14nPxkO23ESEin4tIm9Z0LkRkkYikVLwPra7PLyJP2q6lh0Vktr3jN5tk4eD0IS1ZCfCIMWYgMBa41/b5nwDWGGOCgTW25dbgASCqwnJrPQ//BX40xoQCw7Gek1Z3LkSkO/B7INwYMwTroJpraV3n4gNgTpV1NX5+27XjWmCwbZ83bdfYWjWbZIFj04e0WMaYRGPMbtu/c7BeFLpjPQcf2jb7ELiiUQJsQCISBFwCvFthdWs8Dx2AScB7AMaYImNMJq3wXNi4AW1FxA1oh/WerVZzLowx64HTVVbX9vnnA18YYwqNMcexjkgdXdfxm1OyqG1qkFZHRHoDFwHbgICye1Ns/+3SiKE1lFeAxwBLhXWt8Tz0BVKB921Ncu+KiBet8FwYY04B/wZOYJ0uKMsYs4pWeC6qqO3zn/P1tDklC4emBmnpRKQ98C3woDEmu7HjaWgicimQYozZ1dixNAFuQBjwljHmIiCPlt3MUitbW/x8oA/QDfASkRsaN6om7Zyvp80pWbT6qUFExB1rovjUGLPYtjq5bKZe239TGiu+BjIBuFxEYrE2RU4TkU9ofecBrL+JeGPMNtvyN1iTR2s8FzOA48aYVGNMMbAYGE/rPBcV1fb5z/l62pyShSPTh7RYIiJY26ajjDEvV3hpKXCz7d83A0saOraGZIx50hgTZIzpjfU78LMx5gZa2XkAMMYkASdFpGw20enAQVrhucDa/DRWRNrZfivTsfbrtcZzUVFtn38pcK2IeIpIH6zPFNpe14Ga1R3cIjIPa3t12fQhf2vciBqOiFwMbAAOcLat/o9Y+y2+Anpi/cFcbYyp2snVIonIFOAPxphLRaQTrfA8iMgIrB39HkAM1ilzXGid5+JZ4BqsIwf3AHcA7Wkl50JEPgemYJ2WPRn4M/A9tXx+EXkKuA3r+XrQGLOizuM3p2ShlFKqcTSnZiillFKNRJOFUkopuzRZKKWUskuThVJKKbs0WSillLJLk4VSdthmdr3H9u9uIvJNY8ekVEPTobNK2WGbi2uZbTZTpVolpz2DW6kW5AWgn4jsBY4CA40xQ0TkFqyzeLoCQ4CXsN4cdyNQCMwzxpwWkX5Yp9f3x/qs+d8aYw419IdQ6nxoM5RS9j0BHDPGjAAerfLaEOA6rNM7/w3It03qtwW4ybbNQuB+Y8xI4A/Amw0RtFIXktYslDo/a23PF8kRkSzgB9v6A8Aw2yzB44GvrVMWAeDZ8GEqdX40WSh1fgor/NtSYdmC9fflAmTaaiVKNVvaDKWUfTmAd312tD1z5LiIXA3W2YNFZPiFDE6phqDJQik7jDHpwCYRiQD+VY9DXA/cLiL7gEha0eOAVcuhQ2eVUkrZpTULpZRSdmmyUEopZZcmC6WUUnZpslBKKWWXJgullFJ2abJQSilllyYLpZRSdv0/hTGwiwlQMhoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.xlabel(\"time\")\n",
    "plt.ylabel(\"order parameter\")\n",
    "plt.xlim(t0, t1)\n",
    "plt.ylim(0, 1)\n",
    "plt.plot(sol.ts, sol.observables)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "7d3977cd516aab4e7ac34ec0977b1608119a96f0bf9dd8c065a30d2af58323e0"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}