diff --git a/Gemfile.lock b/Gemfile.lock
index 553ea295..ce747b1e 100644
--- a/Gemfile.lock
+++ b/Gemfile.lock
@@ -12,6 +12,7 @@ GEM
     eventmachine (1.2.7)
     ffi (1.17.0-arm64-darwin)
     ffi (1.17.0-x86_64-darwin)
+    ffi (1.17.0-x86_64-linux-gnu)
     forwardable-extended (2.6.0)
     google-protobuf (4.28.3-arm64-darwin)
       bigdecimal
@@ -19,6 +20,9 @@ GEM
     google-protobuf (4.28.3-x86_64-darwin)
       bigdecimal
       rake (>= 13)
+    google-protobuf (4.28.3-x86_64-linux)
+      bigdecimal
+      rake (>= 13)
     http_parser.rb (0.8.0)
     i18n (1.14.6)
       concurrent-ruby (~> 1.0)
@@ -73,6 +77,8 @@ GEM
       google-protobuf (~> 4.28)
     sass-embedded (1.80.5-x86_64-darwin)
       google-protobuf (~> 4.28)
+    sass-embedded (1.80.5-x86_64-linux-gnu)
+      google-protobuf (~> 4.28)
     terminal-table (3.0.2)
       unicode-display_width (>= 1.1.1, < 3)
     unicode-display_width (2.6.0)
@@ -84,6 +90,7 @@ PLATFORMS
   arm64-darwin-24
   x86_64-darwin-21
   x86_64-darwin-22
+  x86_64-linux-gnu
 
 DEPENDENCIES
   jekyll (~> 4.3.3)
diff --git a/_layouts/sim.html b/_layouts/sim.html
index 56accbdf..bc212b37 100644
--- a/_layouts/sim.html
+++ b/_layouts/sim.html
@@ -82,6 +82,13 @@
         }
       }
     </script>
+    <script type="importmap">
+      {
+        "imports": {
+          "@google/genai": "https://esm.run/@google/genai"
+        }
+      }
+      </script>
     <script async type="text/javascript">
       // If the search string specifies a minified link, query the endpoint for the full link.
       if (searchParams.has("mini")) {
@@ -468,21 +475,20 @@
     </div>
     <div class="message big-popup" id="llm-message">
       <div class="centre_symbol"><i class="fa-solid fa-robot"></i></div>
-      <p>
-        You've clicked the LLM support button! This is an experimental feature
-        that will connect you to a large language model for help with your
-        simulation.
-      </p>
-      <button class="btn-primary message_ok highlight" id="llm_click_open">Open LLM
-        <svg style="color:" xmlns="http://www.w3.org/2000/svg" viewBox="0 0 24 24" width="16" height="16" aria-hidden="true">
-          <path fill="currentColor" d="M14 3h7v7h-2V6.41l-9.29 9.3-1.42-1.42 9.3-9.29H14V3z"></path>
-          <path fill="currentColor" d="M5 5h6v2H7v10h10v-4h2v6H5z"></path>
-        </svg>
-      </button>
+      <p>Equation Generator</p>
+      
+      <div id="ai-link-generator">
+        <textarea id="ai-prompt-input" 
+                  placeholder="Describe a simulation (e.g. 'Slow moving blue waves with high diffusion')"></textarea>
+        
+        <button class="btn-primary highlight" id="ai-generate-btn">
+          Generate
+        </button>
+
+        <div id="ai-status"></div>
+      </div>
+
       <button class="btn-primary message_dismiss" id="llm_click_close">Close</button>
-      <p style="font-size: 0.8em">
-        If you don't see a VisualPDE-tailored chat, please sign in to Google, then open the LLM again.
-      </p>
     </div>
 
     <div id="comboBCs_ui" class="ui hidden">
diff --git a/sim/css/main.css b/sim/css/main.css
index 21bd3f55..c7c3621a 100644
--- a/sim/css/main.css
+++ b/sim/css/main.css
@@ -2007,3 +2007,50 @@ button {
 .colourbarMain {
   filter: brightness(90%);
 }
+
+
+#ai-link-generator {
+  margin-top: 15px;
+  border-top: 1px solid var(--ui-button-border-color);
+  padding-top: 15px;
+  text-align: left;
+}
+
+#ai-prompt-input {
+  width: 100%;
+  height: 80px;
+  padding: 10px;
+  border-radius: 5px;
+  border: 1px solid var(--ui-button-border-color);
+  background-color: var(--input-background-color);
+  color: var(--text-color);
+  font-family: inherit;
+  font-size: 0.9rem;
+  box-sizing: border-box;
+  margin-bottom: 10px;
+  resize: none;
+}
+
+#ai-prompt-input:focus {
+  outline: 2px solid var(--link-color);
+  border-color: transparent;
+}
+
+#ai-generate-btn {
+  width: 100%;
+  margin-bottom: 10px;
+  cursor: pointer;
+}
+
+#ai-status {
+  font-size: 0.8rem;
+  min-height: 1.2rem;
+  color: var(--soft-text-color);
+  margin-bottom: 10px;
+  font-style: italic;
+}
+
+#ai-status a {
+  font-weight: bold;
+  text-decoration: underline;
+}
\ No newline at end of file
diff --git a/sim/llm-context/FAQ.md b/sim/llm-context/FAQ.md
new file mode 100644
index 00000000..6c28b905
--- /dev/null
+++ b/sim/llm-context/FAQ.md
@@ -0,0 +1,99 @@
+# FAQs
+
+## Sharing and citing
+
+### How do I cite VisualPDE in my article?
+The VisualPDE paper is out! You can find it [online](https://doi.org/10.1007/s11538-023-01218-4) for free and cite it as 
+
+Walker, B.J., Townsend, A.K., Chudasama, A.K. et al. VisualPDE: Rapid interactive simulations of partial differential equations. Bulletin of Mathematical Biology 85, 113 (2023).
+
+### Am I allowed to make and share content using VisualPDE.com?
+Yes! We would love to see you making and sharing content using VisualPDE. We ask that any content you create somehow points back to [VisualPDE.com](https://visualpde.com) – a great way to do this is to share a link directly to your simulation! Feel free to use VisualPDE in articles, images, videos, scientific works, social media posts, etc.
+
+### I want to share my cool simulation – how do I do this?
+VisualPDE is designed for sharing. Whether you've fallen in love with one of our examples or you've crafted your own using the VisualPDE interface, you can share the current simulation using a link. Simply open the share sheet by clicking Share (📤) and click 'Copy link' to put a URL on your clipboard that you can share anywhere by pasting! It's that easy!
+
+These links used to be quite long, but now we try to give you a shortened link that works a bit like [tinyURL.com](https://tinyurl.com/app). If you see the 'mini' tag appear next to 'Copy link', you're going to get one of our mini links! Perfect for sharing!
+
+### How do I take a pretty screenshot or video without all the user interface showing up?
+Screenshots and videos are simple in VisualPDE. Simply open the share sheet by clicking Share (📤) and click 'Save image' to download a screenshot to your device, or 'Record clip' to begin a recording that will download to your device when done (we recommend Chrome or Firefox for this). Recordings can be up to 60s long and will stop before then if you pause the simulation or click {{ layout.stop_recording }}
+
+### How do I make my videos compatible with social media and messaging apps?
+Different services accept different types of video, so compatibility is hard. Even worse, different browsers can only create certain types of video, so universal video compatibility is essentially impossible for us at VisualPDE. Fortunately, there are simple ways to convert video so that it is compatible with your favourite services.
+
+If you want a quick solution, <a href="https://ffmpeg-online.vercel.app/?inputOptions=-i&output=VisualPDERecording.mp4&outputOptions=" target="_blank">ffmpeg-online</a> will do this for you. Just drag and drop your video, ignore the rest of the boxes and hit the run button to generate a compatible video file. Even better: your video never leaves your device.
+
+If you prefer a command line solution, then [``ffmpeg``](https://www.ffmpeg.org) is once again your friend. If you have ``ffmpeg`` installed, running ``ffmpeg -i FILENAME output.mp4`` will almost certainly save a compatible copy of your video (``FILENAME``) as ``output.mp4``. Don't just take our word for it - try it out!
+
+### Can I embed VisualPDE in my own site?
+Yes you can! Clicking Share (📤) and then 'Embed' will put an [iframe](https://developer.mozilla.org/en-US/docs/Web/HTML/Element/iframe) element on your clipboard that points to the current simulation. You can specify the complexity of the user interface that you want to be visible, selecting from the full experience, the minimal interface used in our [Visual Stories](/visual-stories), or no user interface at all. If you do this, we'd love to hear from you at [hello@visualpde.com](mailto:hello@visualpde.com) so we can check out your site and your simulation!
+
+In principle, you can have as many embedded simulations on your site as you want. In practice, having lots of simulations running at once might cause stuttering on some devices. You can mitigate this by pausing simulations (VisualPDE consumes very few resources when paused) or by only showing a limited number of simulations onscreen at any one time (simulations that are out of view are automatically paused by your browser). You can also make your simulations less computationally demanding by reducing the number of timesteps per frame or the domain size.
+
+---
+
+## Errors
+Sometimes, mistakes happen. If something has gone wrong, VisualPDE will try its best to describe the error in a helpful way via a pop-up. We try to display a helpful error message to help you identify the problem. More cryptic messages are listed below.
+
+### Undeclared identifier
+```
+Unknown symbol: 'k'
+```
+
+Here, the quantity 'k' has been used somewhere in one of the many free-text inputs in VisualPDE, but it hasn't been defined. The most common cause of this is using a parameter in the **Definitions** tab without defining it in the **Parameters** tab. Check your definitions and parameters to resolve this. Be careful to define parameters as constant numbers or in terms of other parameters – dependence on space, time, or species is not (yet) supported!
+
+### It looks like things have spiralled off to infinity, or the numerical solver is unstable
+```
+Check the equations, or try reducing the timestep, and restart the simulation
+```
+
+The default timestep size (`$\Delta t = 0.1$`) is often too large, try reducing it by an order of magnitude. Alternatively, change the numerical solver.
+
+### Cyclic variables/parameters detected
+```
+Cyclic variables/parameters detected. Please check the definition(s) of a,b.
+```
+
+This error arises when the definitions of your variables or parameters mean that they depend on each other in a cyclic way that can't be resolved. For example, imagine parameters $a$ and $b$ are defined as
+
+$$\begin{align}a &= b,\\ b &= a+1.\end{align}$$
+
+There is no solution to this system of equations, so check your definitions and try again.
+
+### My error looks nothing like any of these
+Sometimes, errors won't look anything like these examples. If this is the case and the error message doesn't help you in resolving it, please follows the steps outlined [below](#error) and help make VisualPDE as stable as possible!
+
+NOTE: VisualPDE sometimes won't warn you about multiple errors that involve the same error message (we're working on it). So, if you fix an error involving 'k', you might not be warned about later errors involving 'k' in the same session. Reloading the page (making sure to have copied your configuration URL first!) is a good way of getting around this for now.
+
+### The buttons don't work, but I want to keep playing!
+Very rarely, the simulation and the user interface may become unresponsive. If this happens: 
+1. Click Help (❓) to bring up the documentation to try to solve your problem.
+2. Reload the page. Sadly, this won't preserve the configuration beyond that specified in the URL.
+If the issue is persistent and you're confident that you've followed all the steps outlined in this guide, please report your issue as described [below](#error). 
+
+### How do I report an issue?
+Sometimes, something might go wrong (e.g. a part of the user interface is doing something strange, or the simulations are not working on your device). If you encounter an issue that can't be resolved using the [documentation](/user-guide), we'd love to hear from you at [bugs@visualpde.com](mailto:bugs@visualpde.com).
+
+If you can, it will help us a lot if you could click <span class='click_sequence'>Settings (🔧) → **Misc.** → **Debug** → **Copy debug**</span> and paste the contents into your email. This will contain information about your current configuration that will help us resolve the problem.
+
+---
+
+## Beyond VisualPDE.com
+
+### I want to use VisualPDE in my teaching/research/outreach – what do I do?
+VisualPDE is designed for customising and sharing. If you're new to the site, we recommend that you check out our [creator resources](/create) to get an idea of how to craft your own simulation, and make use of [link sharing](#linkSharing) to permanently save customised simulations.
+
+For instance, a popular way to create custom teaching materials is to simply include links to customised VisualPDE simulations in more traditional teaching materials, similar to the [examples](/explore) on the site.
+
+If you want to do more than this allows, we'd love to hear from you at [hello@visualpde.com](mailto:hello@visualpde.com) so that we can help bring VisualPDE into your teaching, research, or outreach activities.
+
+### I want to use VisualPDE for my business - can I?
+Almost certainly! However, we want to make sure that we do knowledge exchange properly, so do get in touch with us at [hello@visualpde.com](mailto:hello@visualpde.com) so that we can make sure you'll get the best out of VisualPDE.
+
+### I've got a feature request – who do I contact?
+We're always looking to improve VisualPDE and would love to hear from you at [hello@visualpde.com](mailto:hello@visualpde.com) if you think we can improve something about the experience. If you'd like to suggest a brand new PDE for the website, we're looking for examples that are qualitatively different to what is already on the site. That being said, if you've found a cool parameter set for a system that we haven't noted already, do let us know!
+
+If you can't seem to cast your PDE in a form that VisualPDE can solve, we encourage you to look at the [library of examples](/explore) to see how various types of PDE can be thrown into VisualPDE, which might provide some inspiration. If you still can't get VisualPDE to solve your system but you think it should/could be possible, we'd love to hear from you.
+
+### Can I modify the source code of VisualPDE and host my own version?
+In short: yes! We've made VisualPDE open source for a reason and want to see as many people using it as possible. If you want to ship your own version of the source code, we ask that you follow the licences found in our [repository](https://github.com/Pecnut/visual-pde) and provide due credit to [VisualPDE.com](https://visualpde.com). If you're in any doubt about your specific case, do send us an email at [hello@visualpde.com](mailto:hello@visualpde.com) and we'll happily chat with you! We're all about making VisualPDE accessible to as many people as possible.
\ No newline at end of file
diff --git a/sim/llm-context/advanced-options.md b/sim/llm-context/advanced-options.md
new file mode 100644
index 00000000..1ae86d0e
--- /dev/null
+++ b/sim/llm-context/advanced-options.md
@@ -0,0 +1,681 @@
+# Complete documentation
+
+Almost everything in VisualPDE is customisable. Here, we describe the basic functionality of every option that can be found in the menus of VisualPDE.
+
+## Equations ($f(x)$)
+
+VisualPDE is all about solving equations. In the Equations pane, you can view and define the problem that VisualPDE will solve for you in your browser, complete with initial and boundary conditions. More advanced settings, including variable renaming, can be found under [**Settings**](#settings).
+
+### Edit
+
+Customise all the terms in the PDEs that you would like to solve using natural syntax. See our discussion of [valid expressions](#writing-valid-expressions) for helpful examples that will guide you in posing your own PDE system. Typing in any of the fields will highlight the corresponding term in the typeset PDE above.
+
+- #### Typeset
+
+  Have VisualPDE typeset the specified equations, making use of all the defined diffusion coefficients, functions and parameters. Terms will not be substituted in if they are constants that are not 0 or 1. Toggle this off to see the format of the equations that VisualPDE can interpret. 
+  JSON key: `typesetCustomEqs`, boolean. Default: `true`
+
+- #### $D_u$, $D_v$, $D_w$, ...
+
+  Set the diffusion coefficients of all the species in the simulation. When **Cross diffusion** is enabled, you can also set interaction terms, which are written $D_{uv}$ etc. These can be functions of space ($x$, $y$), time ($t$), any of the unknowns ($u$, $v$, $w$, $q$), the size of the domain ($L$, $L_x$, $L_y$), the images ($I_S$, $I_T$) and any quantities defined in **Parameters**. See our discussion of [valid expressions](#writing-valid-expressions) for valid syntax and a list of available in-built functions.
+
+  Advanced users seeking diagonal anisotropic diffusion tensors (i.e. distinct diffusion coefficients in the coordinate directions) can define two coefficients at once by separating definitions with a semicolon. For example, `1;2` sets diffusion coefficients of `1` and `2` in the $x$ and $y$ directions, respectively.
+  JSON key: `diffusionStr_i_j`: definitions, string. Defaults: `diffusionStr_1_1`: `"1"`, `diffusionStr_2_2`: `"2"`, `diffusionStr_3_3`: `"0"`.
+
+- #### $f_u$, $f_v$, $f_w$, ...
+
+  Define the inhomogeneities in the equations. These can be functions of space ($x$, $y$), time ($t$), any of the unknowns ($u$, $v$, $w$, $q$), the size of the domain ($L$, $L_x$, $L_y$), the images ($I_S$, $I_T$), and any quantities defined in **Parameters**. See our discussion of [valid expressions](#writing-valid-expressions) for valid syntax and a list of available in-built functions.
+
+  Advanced users can also make careful use of `RAND`, a uniformly random value in $[0,1]$, and `RANDN`, a normally distributed random number with unit variance and zero mean. This converts the equations into [stochastic partial differential equations](https://en.wikipedia.org/wiki/Stochastic_partial_differential_equation), which should only be solved using the Forward Euler timestepping scheme. Both `RAND` and `RANDN` require manually dividing by `sqrt(dt)` in non-algebraic equations so that the scheme resembles the [Euler-Maruyama method](https://en.wikipedia.org/wiki/Euler–Maruyama_method). The solution under other timestepping schemes is undefined.
+
+  For convenience, we define `WhiteNoise` to be a normally distributed random variable with unit variance and zero mean, scaled by $1/\sqrt{\Delta t\Delta x^N}$ where $N=1,2$ is the spatial dimension. This scales appropriately with timestep and spatial step and is suitable for direct use in PDEs, such as in our [stochastic example](/nonlinear-physics/stochastic-pdes). You can use up to 4 independent WhiteNoise terms using the syntax `WhiteNoise_1`, `WhiteNoise_2`, `WhiteNoise_3` and `WhiteNoise_4`.
+
+  JSON key: `reactionStr_1` ... `reactionStr_4`: definitions, string. Defaults: `reactionStr_1`: `"u^2*v - (a+b)*u"`, `reactionStr_2`: `"-u^2*v + a*(1 - v)"`, `reactionStr_3`: `"0"`.
+
+### Parameters
+
+This menu contains a list of all the user-specified values that can be used throughout VisualPDE. New parameters can be defined using the empty input field at the bottom of the list of parameters. Parameters can depend on one another, but their definitions cannot be cyclic.
+JSON key: `kineticParams`, definitions, string. Default: `"a=0.037 in [0,0.1];b=0.06 in [0.04,0.1];",`
+
+#### Basics
+
+The basic syntax for defining a parameter is
+
+```
+name = value
+```
+
+which will make the quantity `name` available to the simulation. You can then freely change `value`, which will instantly propagate throughout VisualPDE. If you try to use a name that clashes with an internal variable (some of which are only found under the hood of VisualPDE), a warning will appear to inform you of this. Parameters can be removed by deleting the text that defines them. You can even choose a `name` that includes subscripts, such as `k_1u`. This will be interpreted as $k_{1u}$ automatically by VisualPDE.
+
+#### Sliders
+
+The more advanced syntax
+
+```
+name = value in [start,step,stop]
+```
+
+creates a slider for your variable, ranging between the `start` and `stop` values in increments of `step`. The `step` parameter can be omitted and VisualPDE will choose a step automatically. For example,
+
+```
+a = 0.5 in [0,1]
+```
+
+creates a slider that ranges between 0 and 1, with initial value 0.5 and an automatically determined step size. Parameters with sliders cannot be defined in terms of other parameters.
+
+The configuration of a slider (value, start, step, stop) can be updated by modifying the relevant parts of the expression that defines it. Sliders can be removed by deleting `in ...` from the parameter definition, and will be removed automatically when the associated parameter is removed.
+
+### Boundary conditions
+
+Boundary conditions can be specified for any species in the simulation. The following boundary conditions are available:
+
+- Periodic
+- [Dirichlet](https://en.wikipedia.org/wiki/Dirichlet_boundary_condition) (e.g. $u|_{\partial\Omega} = 0$)
+- [Neumann](https://en.wikipedia.org/wiki/Neumann_boundary_condition) (e.g. $\frac{du}{dn}|_{\partial\Omega} = 0$)
+- [Robin](https://en.wikipedia.org/wiki/Robin_boundary_condition) (e.g. $(u + \frac{du}{dn})|_{\partial\Omega} = 0$)
+
+Boundary conditions that allow you to specify values can be functions of space ($x$, $y$), time ($t$), any of the unknowns ($u$, $v$, $w$, $q$), the size of the domain ($L$, $L_x$, $L_y$), the images ($I_S$, $I_T$) and any quantities defined in **Parameters**. Robin boundary conditions are the only type supported that allow you to use an unknown in the specification of its own boundary condition. See our discussion of [valid expressions](#writing-valid-expressions) for valid syntax and a list of available in-built functions.
+
+An additional option, **Mixed...**, is also available, which allows you to specify different types of boundary condition on the Left, Right, Top and Bottom sides of rectangular domains.
+
+The easiest way to do this is using the graphical interface by clicking 🧊 beside the boundary condition selector (unavailable during surface plots). Click the various boundaries to set the desired condition on each. These conditions can also be specified directly as a string, e.g.
+
+```
+Left: Dirichlet = 0; Right: Neumann = 1; Top: Robin = u; Bottom: Dirichlet = sin(x)
+```
+
+for the species $u$ would specify $u = 0$ on the left boundary, $\frac{du}{dn} = 1$ on the right boundary, $\frac{du}{dn} = u$ on the top boundary and $u = \sin(x)$ on the bottom boundary. Sides can be specified in any order and are case sensitive. Omitting any side will default to periodic boundary conditions (beware, this may have unexpected results if the matching side is not also periodic; using the graphical interface prevents this).
+
+An additional type of condition, 'Ghost', can also be specified with Mixed boundary conditions. This advanced option pushes VisualPDE to its limits, overriding the value of the [ghost nodes](https://kyleniemeyer.github.io/ME373-book/content/bvps/finite-difference.html#using-central-differences-for-derivative-bcs) used in the spatial discretisation of the PDE, and should be used with caution. We make use of this option in our Visual Story on [virus transmission](/visual-stories/airborne-infections) to effectively double the size of the computational domain in one direction. This must be toggled on in <span class='click_sequence'>Settings (🔧) → **Misc.**</span>
+JSON key: `boundaryConditions_i`, one of (`periodic`, `dirichlet`, `neumann`, `robin`). Defaults: `boundaryConditions_1`: `periodic`, `boundaryConditions_2`: `periodic`.
+JSON key: `dirichletStr_i`, definition, string.
+JSON key: `neumannStr_i`, definition, string.
+JSON key: `robinStr_i`, definition, string.
+
+### Initial conditions
+
+Initial conditions can be specified for any species in the simulation. They can be functions of space ($x$, $y$), the size of the domain ($L$, $L_x$, $L_y$), the images ($I_S$, $I_T$), the random quantity `RAND`, a uniformly random value in $[0,1]$, the random quantity `RANDN`, a normally-distributed random number with unit variance and zero mean, and any quantities defined in **Parameters**. See our discussion of [valid expressions](#writing-valid-expressions) for valid syntax and a list of available in-built functions.
+JSON key: `initCond_i`, definition, string. Defaults: `initCond_1`: `"0"`, `initCond_2`: `"1"`.
+
+### Advanced options
+
+Configure additional equation-related settings, including the number of species and the type of terms that will be included.
+
+- #### \# species
+
+  Specify the number of unknowns (1, 2, 3, or 4) in the simulation.
+  JSON key: `numSpecies`, one of `1`, `2`, `3` or `4`. Default: `2`.
+
+- #### \# algebraic
+
+  Choose how many equations you want to be in algebraic form in systems with cross diffusion enabled. The equations will be put in algebraic form in reverse order, e.g. a 4-species system with 1 algebraic species will convert the final equation to be algebraic.
+  JSON key: `numAlgebraicSpecies`, one of `0`, `1`, `2`, or `3`. Default: `0`.
+
+- #### Species names
+
+  Specify custom names for the species in VisualPDE. When changing "# species" to a higher number previously set, new species will be named by default to `SPECIES2`, `SPECIES3` and `SPECIES4`, for 2, 3, 4 species, respectively. Names can be multi-character and can include letters, numbers, and underscores, but must each be a single 'word'. For example, `T_01` is a valid name (rendered as $T_{01}$) whilst `T 01` is not. Space or commas can be used to separate names in the list. Certain names are reserved under the hood, such as `H` for the Heaviside function, but VisualPDE will warn you if you attempt to use a reserved name. VisualPDE will automatically substitute the names of old species everywhere in the simulation and interface.
+  JSON key: `speciesNames`, space separated string. Default: `"u v"`
+
+- #### Cross diffusion
+
+  Enable cross diffusion in systems with 2 or more species, enabling simulation of a wide range of systems in which a species can depend on the gradient of another.
+  JSON key: `crossDiffusion`, boolean. Default: `false`.
+
+- #### Scales
+
+  Set per-equation timescales (multiplying any time derivatives) $\tau_u$, $\tau_v$, $\tau_w$, $\tau_q$ to enable simpler entry of some types of systems. For algebraic equations, these quantities are no longer timescales, but retain their notation and function as per-equation scale factors. They can be functions of space ($x$, $y$), time ($t$), any of the unknowns ($u$, $v$, $w$, $q$), their gradients ($u_x$, $u_y$, etc.), the size of the domain ($L$, $L_x$, $L_y$), the images ($I_S$, $I_T$), and any quantities defined in **Parameters**. See our discussion of [valid expressions](#writing-valid-expressions) for valid syntax and a list of available in-built functions.
+
+  Importantly, **timescales must be non-zero**. Setting timescales to zero will result in singularities and are equivalent to large diffusion coefficients, large timesteps, or fast kinetic terms.
+  JSON key: `timescales`, enable scales, boolean. Default: `false`.
+  JSON key: `timescale_i`, definition, string.
+
+---
+
+## Views (📚)
+
+There are often multiple ways to visualise a solution to a PDE. In the Views pane, you can select from and customise a range of example-specific display options, or create your own. Everything you customise will be saved in the current View. If you share your simulation via a link, your Views will be sent along too.
+
+### New (+)
+
+Create a new view with a placeholder name from the current view configuration.
+
+### Rename
+
+Edit the name of the current View, enclosing any mathematics in `$` tags. You can even use emoji.
+
+### Delete
+
+Delete the currently selected View. Only visible if there are at least two views.
+
+### Expression
+
+Choose the expression that you want to be used to colour the domain, which can be any function of the species solved for, as well as space, time, and user-defined parameters. Often, this is either $u$, $v$, $w$ or $q$. Explicitly, this can be a function of space ($x$, $y$), time ($t$), any of the unknowns ($u$, $v$, $w$, $q$) and their gradients ($u_x$, $u_y$, etc.), the size of the domain ($L$, $L_x$, $L_y$) and the images ($I_S$, $I_T$).
+JSON key: `whatToDraw`, `whatToPlot`, definition, string. Defaults: `whatToDraw`: `"u"`, `whatToPlot`: `"u"`.
+
+### Plot type
+
+Choose from three types of plot: **line**, **plane** or **surface**. Any simulation can be viewed as any plot type.
+
+Line plots are the default plot type for 1D domains. Cubic splines are used to interpolate between nodes of the computational domain for smooth plotting. This may lead to transient oscillations appearing near discontinuities in the solution.
+
+Surface plots are constructed by using the chosen **Expression** as a height map, the limits of the colour axis and the **Height scale** parameter.
+JSON key: `plotType`, one of (`line`, `plane`, `surface`).
+
+### Colour
+
+Customise everything about the colours used to display the solution.
+
+- #### Colour map
+
+  Set the current colour map being used to convert **Expression** into a colour value. Use the dropdown to select from the available options. We have tried to cater for everyone in these options but, if you find that no colour map is available that allows you to easily distinguish between values, please let us know at [hello@visualpde.com](mailto:hello@visualpde.com) so that we can add a more appropriate map.
+  JSON key: `colourmap`, one of (`redGreen`, `fireOnTerrain`, `splitscreenFires`, `cyclic`, `pride`,
+  `terrain`, `squirrels`, `chemicalBlue`, `chemicalYellow`, `chemicalGreen`,
+  `spooky`, `retro`, `greyscale`, `urbanFlooding`, `BlackGreenYellowRedWhite`,
+  `viridis`, `turbo`, `blue-magenta`, `diverging`, `thermal`,
+  `snowghost`, `midnight`, `lavaflow`, `ice`, `pastels`,
+  `foliage`, `water`, `blue`). Default: `viridis`
+
+- #### Min/Max value
+
+  Set the limits of the colour map that transforms **Expression** into colour. If viewing a surface plot, this also impacts the height of the surface. These can depend on any quantities defined in **Parameters**.
+  JSON key: `minColourValue`, string. Default: `"0"`.
+  JSON key: `maxColourValue`, string. Default: `"0.37"`.
+
+- #### Reverse
+
+  Reverse the direction of the current colour map.
+  JSON key: `flippedColourmap`, boolean.
+
+- #### Colour bar
+
+  Toggle the display of the current colour bar and limits.
+  JSON key: `colourbar`, boolean.
+
+- #### Snap range
+
+  Click to instantly snap **Min value** and **Max value** to the current minimum and maximum of **Expression** in the domain.
+
+- #### Auto snap
+  Toggle the automatic snapping of the colour map limits. This can be very useful if you don't know the range in which **Expression** will fall, especially if it is changing frequently.
+  JSON key: `autoSetColourRange`, boolean.
+
+### Contours
+
+Toggle the rendering of contours on top of the simulation display. The number, colour, and sensitivity of these (equally spaced) contours can be configured in the **Contours** menu that appears when contours are enabled.
+JSON key: `contours`, boolean.
+
+- ##### Colour
+
+  Specify the colour in which contours will be rendered.
+  JSON key: `contourColour`, hex string.
+
+- ##### Number
+
+  Specify the number of contour lines to plot. These are drawn at equally spaced intervals between the minimum and maximum of the current colourmap e.g. specifying one contour will draw a line at the midpoint of the colourmap.
+  JSON key: `contourNum`, integer.
+
+- ##### Threshold
+  Set the relative numerical threshold, between 0 and 1, within which a contour will be detected. Smaller values result in more precise, thinner contours. Larger values may be needed to detect contours in solutions that vary rapidly in space. Under the hood, VisualPDE checks if a given output pixel is within this threshold of a given contour value, with all values normalised by the range of the colourbar. If the range of the colourbar is 0, the behaviour of contours is undefined.
+  JSON key: `contourEpsilon`, float.
+
+### Lighting
+
+Toggle lighting effects, which adds reflections and shadows to the solution. This often adds a fluid-like character to a simulation. We make use of the [Phong reflection model](https://en.wikipedia.org/wiki/Phong_reflection_model). Details of the filter, including its strength and the orientation of the simulated light, can be specified in the **Lighting** menu that appears when lighting is enabled. Some lighting effects may appear slightly pixellated on some devices (typically Android tablets and iPadOS devices), though increasing the grid refinement will mitigate this.
+JSON key: `emboss`, boolean.
+
+- #### Smoothness
+
+  Configure the simulated smoothness of the surface. Low values will result in sharp shadows and reflections, whilst larger values will produce an apparently smoother surface.
+  JSON key: `embossSmoothness`, float.
+
+- #### Ambient
+
+  Set the level of background light, which does not contribute to reflections or shadows.
+  JSON key: `embossAmbient`, float.
+
+- #### Diffuse
+
+  Set how much the light source impacts on the brightness of the surface.
+  JSON key: `embossDiffuse`, float.
+
+- #### Specular
+
+  Specify the intensity of specular (shiny) reflections.
+  JSON key: `embossSpecular`, float.
+
+- #### Precision
+
+  Set the precision (also known as the specular exponent) of the specular reflections. Low values widely scatter light, whilst large values result in precise highlights.
+  JSON key: `embossShiny`, float.
+
+- #### Inclination
+
+  Set the angle of inclination of the light source. Setting this to 0 places the light source directly overhead.
+  JSON key: `embossTheta`, float.
+
+- #### Direction
+  Specify the in-plane direction of the light source, rotating any shadows and highlights.
+  JSON key: `embossPhi`, float.
+
+### Overlay
+
+Toggle the display of an overlay. The expression, colour, and threshold used in displaying the overlay can be specified in the **Overlay** menu that appears when the overlay is enabled.
+JSON key: `overlay`, boolean.
+
+- #### Colour
+
+  Specify the colour in which overlays will be rendered.
+  JSON key: `overlayColour`, hex string.
+
+- #### Expression
+
+  Set an expression whose zero set defines a curve to be displayed in the domain. This can be a function of space ($x$, $y$), time ($t$), any of the unknowns ($u$, $v$, $w$, $q$), their gradients ($u_x$, $u_y$, etc.), the size of the domain ($L$, $L_x$, $L_y$), the images ($I_S$, $I_T$), and any quantities defined in **Parameters**. See our discussion of [valid expressions](#writing-valid-expressions) for valid syntax and a list of available in-built functions.
+  JSON key: `overlayExpr`, definition, string.
+
+- #### Threshold
+  Set the relative numerical threshold, between 0 and 1, within which the zero set will be detected. Smaller values result in more precise, thinner curves. Larger values may be needed to detect curves in large spatial domains or for expressions that vary rapidly in space. Under the hood, VisualPDE checks if a given pixel is in the zero set to within this threshold.
+  JSON key: `overlayEpsilon`, float.
+
+### 3D options
+
+When viewing surface plots, this menu will appear to allow you to customise aspects of the display.
+
+- #### Custom surface
+
+  Toggle the rendering of the solution on a custom, user-defined surface. The surface $z(x,y)$ is specified in **Surface $z$**, which appears when a custom surface is enabled. This definition can be a function of space ($x$, $y$), time ($t$), any user-defined parameters, any of the unknowns ($u$, $v$, $w$, $q$) and their first derivatives, the size of the domain ($L$, $L_x$, $L_y$) and the images ($I_S$, $I_T$).
+  JSON key: `customSurface`, boolean.
+
+- #### Height scale
+
+  The scale factor used when plotting surfaces or lines. Changing this parameter effectively makes the vertical variation appear more/less prominent, though it does not alter the plotted colours. Must be a numerical value.
+  JSON key: `threeDHeightScale`, float.
+
+- #### View $\theta$/$\phi$
+
+  [Euler angles](https://en.wikipedia.org/wiki/Euler_angles) specifying the current 3D viewpoint, with $\theta\in[0,\pi]$ and $\phi\in[0,2\pi]$. You can manipulate these values either by inputting new values, or see them update as you rotate the viewpoint with your pointer (click and drag). As Euler angles [don't do a perfect job](https://en.wikipedia.org/wiki/Gimbal_lock) of describing orientations, you may (rarely) find that a viewpoint loaded in from a URL isn't quite what you expected.
+  JSON key: `cameraTheta`, float.
+  JSON key: `cameraPhi`, float.
+
+- #### Zoom
+  The current zoom level of the 3D view. You can manipulate the value directly or, on touch devices, by 'pinching to zoom'.
+  JSON key: `cameraZoom`, float.
+
+### Line options
+
+When viewing line plots, this menu will appear to allow you to customise aspects of the display.
+
+- #### Max height
+
+  The maximum height of the plotted line, measured in units of $L$. Changing this parameter effectively makes the vertical variation appear more/less prominent. Must be a numerical value.
+
+- #### Thickness
+  The thickness of the plotted line relative to the default. Must be a numerical value.
+
+### Vector field
+
+Toggle the rendering of a vector field on top of the simulation. The definition, colour, density, and size of the vectors can be fully customised. This option is often used to visualise flows or fluxes.
+JSON key: `vectorField`, boolean.
+
+- #### Colour
+
+  Specify the colour in which arrows will be rendered.
+  JSON key: `arrowColour`, hex string.
+
+- #### $x$, $y$ component
+
+  Specify the $x$ and $y$ components of the vector field. These components can be functions of space ($x$, $y$), time ($t$), any user-defined parameters, any of the unknowns ($u$, $v$, $w$, $q$) and their first derivatives, the size of the domain ($L$, $L_x$, $L_y$) and the images ($I_S$, $I_T$).
+  JSON key: `arrowX`, definition, string.
+  JSON key: `arrowY`, definition, string.
+
+- #### Density
+
+  Specify the density of the rendered arrows, normalised between 0 and 1.
+  JSON key: `arrowDensity`.
+
+- #### Scaling
+
+  Configure how the rendered arrows will be scaled. **_Relative_** normalises arrows relative to a maximum length that you can specify. **_Auto_** normalises arrows relative to the largest computed arrow at the current instant. **_None_** renders all arrows at the same scale.
+  JSON key: `arrowScale`, one of (`relative`, `auto`, `none`).
+
+- #### Max length
+  Specify the constant length by which all arrows will be normalised. Must be a mathematical expression that is not written in terms of any parameters or user-defined quantities.
+  JSON key: `arrowLengthMax`, definition, string.
+
+### Time series
+
+Toggle the display of a live time series graph of the simulation. This can show either the value of the **Expression** at a single (configurable) position, or the integral of the expression over the domain.
+JSON key: `probing`, boolean.
+
+- #### Type
+
+  Choose whether to plot a single-point sample or an integral of the configured **Expression**.
+  JSON key: `probeType`, one of (`point` or `integral`).
+
+- #### Expression
+
+  Specify the expression to be sampled or integrated. This can be a function of space ($x$, $y$), time ($t$), any user-defined parameters, any of the unknowns ($u$, $v$, $w$, $q$) and their first derivatives, the size of the domain ($L$, $L_x$, $L_y$) and the images ($I_S$, $I_T$). Boundary conditions may not be accurately reflected in computed values.
+  JSON key: `probeFun`, definition, string. Default: `"u"`.
+
+- #### $x$, $y$
+
+  Specify the position at which to sample the **Expression**. This can be a function of time ($t$), any user-defined parameters, and the size of the domain ($L$, $L_x$, $L_y$). This can be automatically filled in using the ⌖ button.
+  JSON key: `probeX`, definition, string.
+  JSON key: `probeY`, definition, string.
+
+- #### Duration
+
+  Specify the length of the time series (in units of time). The series will be updated at regular intervals. For the best experience, this should be around 10-400 times the timestep times the number of timesteps per frame.
+  JSON key: `probeLength`, float. Default: `3000`.
+
+---
+
+## Settings (🔧)
+
+Here you can edit a wide range of settings, from the size of the brush to the timestep of the simulation.
+
+### Brush
+
+VisualPDE allows you to interact directly with simulations via a brush by simply clicking/pressing on the domain. The brush paints values onto the discrete representation of the domain, which act like initial conditions for the rest of the simulation. Using the alternate mouse button (often the right mouse button) will negate the value of the brush. The brush can be used to create complex initial conditions, perturbations, or even to (crudely) simulate the effect of a moving boundary.
+
+- #### Enable brush
+
+  Enable or disable the brush. Most simulations will have the brush enabled by default.
+  JSON key: `brushEnabled`, boolean.
+
+- #### Type
+
+  Change the shape of the brush, choosing between **Disk**, **Square**, **Horizontal line** and **Vertical line**. An additional option, **Custom**, allows you to define a custom shape in the **_Indicator_** field by typing in an expression. The brush will draw wherever the expression is positive. Expressions can be a function of space ($x$, $y$), the brush coordinates ($xB$, $yB$), time ($t$), any user-defined parameters, any of the unknowns ($u$, $v$, $w$, $q$), the size of the domain ($L$, $L_x$, $L_y$), and the images ($I_S$, $I_T$).
+  JSON key: `brushType`, one of (`circle`, `square`, `hline`, `vline`, `custom`).
+
+- #### Value
+
+  Change the **value** that you are painting. This can be a function of space ($x$, $y$), time ($t$), any user-defined parameters, any of the unknowns ($u$, $v$, $w$, $q$), the size of the domain ($L$, $L_x$, $L_y$), the images ($I_S$, $I_T$), `RAND`, a uniformly random value in $[0,1]$, and `RANDN`, a normally-distributed random number with unit variance and zero mean. Using the alternate mouse button (often the right mouse button) will negate the value of the brush (unavailable on touch devices).
+  JSON key: `brushValue`, definition, string.
+  JSON key: `brushAction`, one of (`replace`, `add`).
+
+- #### Radius
+
+  Change the brush size, measured on the same scale as the domain size. This can even be a function of space ($x$, $y$), time ($t$), any user-defined parameters, any of the unknowns ($u$, $v$, $w$, $q$), the size of the domain ($L$, $L_x$, $L_y$) and the images ($I_S$, $I_T$).
+  JSON key: `brushRadius`, definition, string.
+
+- #### Species
+  Set the **species** ($u$, $v$, $w$, $q$) you are painting.
+
+### Domain
+
+- #### Dimension
+
+  Choose between a 1D or a 2D computational domain. Switching to 1D effectively removes the $y$ dimension from the simulation. Make sure that any expressions you've defined don't contain a $y$ after moving to 1D.
+  JSON key: `dimension`, one of (`"1"` or `"2"`).
+
+- #### Largest side
+
+  Change the largest side $L$ of the domain. Must be a mathematical expression that is not written in terms of any other parameters or user-defined quantities. If this is not an integer multiple of the space step $\Delta x=\Delta y$, the domain will extend to $\left\lfloor L / \Delta x\right\rfloor\Delta x$ in practice. This guarantees the use of a precise space step.
+  JSON key: `domainScale`, definition, string. Default: `"1000"`.
+
+- #### Space step
+
+  Set the space step $\Delta x=\Delta y$ used in discretising the domain. You may have to decrease the timestep $\Delta t$ in order to maintain numerical stability if you decrease the space step (as discussed in the [timestepping guide](/user-guide/solver#timestepping)). Must be a mathematical expression that is not written in terms of any other parameters or user-defined quantities.
+  JSON key: `spatialStep`, definition, string. Default: `"3"`.
+
+- #### Min. $x$, $y$
+
+  Set the minimum values of $x$ and $y$ in the simulation. This amounts to translating the simulation domain by $(x_{min}, y_{min})$. A common use of this is to centre the simulation domain at $(0,0)$ rather than $(L_x/2, L_y/2)$. Can be a function of $L_x$, $L_y$ and any of the user-defined parameters.
+  JSON key: `minX`, definition, string.
+  JSON key: `minY`, definition, string.
+
+- #### Fill screen
+
+  Choose between the domain filling the display (often recommended) or being forced to be square (vital for guaranteeing the shape of the domain). When selected, the largest edge of the display will correspond to the length set in **Largest side**, so that full-screen simulations can be thought of as cropped, zoomed counterparts to square simulations.
+  JSON key: `squareCanvas`, boolean.
+
+- #### Custom
+
+  Toggle the use of a custom domain $\Omega$ that is determined implicitly from a user-set expression.
+  JSON key: `domainViaIndicatorFun`, boolean.
+
+- #### Ind. fun (indicator function)
+  Define the domain implicitly by setting a boolean (e.g. $x<0.5$) or a simple expression (e.g. $x-0.5$), where (strict) positivity identifies the interior of the domain. This can be a function of space ($x$, $y$), time ($t$), any of the unknowns ($u$, $v$, $w$, $q$), the size of the domain ($L$, $L_x$, $L_y$) and the images ($I_S$, $I_T$). In order to allow VisualPDE to correctly apply boundary conditions, the edge of the computational domain is always counted as being outside of the custom domain.
+  JSON key: `domainIndicatorFun`, definition, string.
+
+### Timestepping
+
+- #### Steps/frame
+
+  Set how many timesteps will be performed every time your browser requests a frame from VisualPDE. This setting effectively allows you to speed up/slow down the simulation without altering the timestep, though large values may cause some stuttering on some devices. Must be a numerical value.
+  JSON key: `numTimestepsPerFrame`, integer.
+
+- #### Timestep
+
+  Set the timestep $\Delta t$ used in the solver. You may have to increase the spatial step $\Delta x$ in order to maintain numerical stability if you increase the timestep (as discussed in the [solver guide](/user-guide/solver)). Must be a numerical value.
+  JSON key: `dt`, float. Default: `0.1`.
+
+- #### Scheme
+
+  Select one of various timestepping schemes. [Forward Euler](https://en.wikipedia.org/wiki/Euler_method) is the fastest but least accurate; the [Midpoint Method](https://en.wikipedia.org/wiki/Midpoint_method) and [Runge-Kutta 4](https://en.wikipedia.org/wiki/Runge–Kutta_methods) improve upon the accuracy and stability of Forward Euler, though are associated with increased computational cost. [Adams-Bashforth 2](https://en.wikipedia.org/wiki/Linear_multistep_method#Two-step_Adams–Bashforth) is more accurate but less stable than Forward Euler. Use of higher accuracy schemes may require a reduction of Steps/frame to reduce stuttering due to increased computational load. When solving stochastic partial differential equations, only Forward Euler is supported.
+  JSON key: `timesteppingScheme`, one of (`Euler`, `Mid`, `RK4`, `AB2`). Default: `Euler`.
+
+- #### Elapsed time
+
+  Show/hide the elapsed time since the simulation was loaded/reset.
+  JSON key: `timeDisplay`, boolean.
+
+- #### Auto pause
+
+  Set the simulation to be automatically paused when the time ($t$) passes a custom threshold (which can depend on user-defined parameters), which can be configured when this option is enabled. You can resume an auto-paused simulation by pressing ▶️.
+  JSON key: `autoPause`, boolean.
+  JSON key: `autoPauseAt`, float.
+
+- #### Performance mode
+  Reduces the quality of the display in order to boost simulation performance. Different simulations and devices will benefit differently from this setting. On average, we see an increase of around 6fps on a Mac Mini M1 on the 'BenchmarkingFast' preset in 4K resolution.
+  JSON key: `performanceMode`, boolean.
+
+### Images
+
+- #### $I_S$, $I_T$
+  Define the scalar fields $I_S$ and $I_T$, which are derived from images that you can upload by clicking on the current image. Via the symbols `I_S` and `I_T` throughout VisualPDE, you can access the average RGB value of each image at each point in space, effectively treating them as greyscale. Advanced users can access the individual RGBA channels via `I_SR`, `I_SG`, etc. VisualPDE will stretch images so that they cover the domain edge-to-edge. Note that this does not respect **Implicit**. The defaults draw from images of [Sofya Kovalevskaya](https://en.wikipedia.org/wiki/Sofya_Kovalevskaya) and [Alan Turing](https://en.wikipedia.org/wiki/Alan_Turing).
+  JSON key: `imagePathOne`, path string.
+  JSON key: `imagePathTwo`, path string.
+
+### Checkpoints
+
+VisualPDE supports checkpoints, which allow you to save the state of a simulation at the touch of a button. This allows you to instantly return to a previous solution state - very handy if you've crafted the perfect initial condition by painting with the brush. Revert to a checkpoint by pressing Restart (🔄)
+
+- #### Enable checkpoints
+
+  Toggle the use of checkpoints. When enabled, resetting the simulation will revert to a saved checkpoint (if one exists) instead of using any initial conditions defined alongside the equations.
+  JSON key: `resetFromCheckpoints`, boolean.
+
+- #### Set
+
+  Save the current simulation as a checkpoint.
+
+- #### Export
+
+  Click to export the last checkpoint as a file to your device, which can be shared and uploaded to the site and used as a checkpoint. If no checkpoint exists, one will be created. By default, the file will be called `VisualPDEState`.
+
+- #### Import
+
+  Import a checkpoint from a VisualPDE file. By default, these are called `VisualPDEState`.
+
+- #### Resize
+  Specify how a checkpoint should be resized to fit the current simulation domain. 'Stretch' will stretch the checkpoint so that it fills the current domain, but will not preserve the aspect ratio in general. "Crop" will crop the checkpoint whilst preserving the aspect ratio, but may result in some information not being used.
+  JSON key: `resizeCheckpoints`, one of (`stretch`, `crop`).
+
+### Misc
+
+- #### Background
+
+  Set the background colour of the simulation window, which you will see often when using **Surface Plot** or **Implicit**. In implicit domains, the exterior of the domain adopts this colour.
+  JSON key: `backgroundColour`, hex string.
+
+- #### Integrate
+
+  Display the integral of **Expression** over the domain. This integral is coarsely approximated by a simple Riemann sum, with accuracy improving with mesh refinement. This can be used to track the numerical evolution of quantities that should be conserved in a system.
+  JSON key: `integrate`, boolean.
+
+- #### Interpolate
+
+  Use this option to force the use of manual, unoptimised filtering in place of device-default interpolation of displayed colours. This toggle is not available on devices that do not support interpolation by default; in this case, manual interpolation is always enabled.
+  JSON key: `forceManualInterpolation`, boolean.
+
+- #### Set seed
+
+  Set the seed of the (pseudo)random number generator used to assign values to `RAND` and `RANDN` in all free-text fields in the VisualPDE interface. Note that `RAND` and `RANDN` always vary in space. A specific random seed in the form of a numerical value can be configured (default 0).
+  JSON key: `setSeed`, boolean.
+  JSON key: `randSeed`, integer.
+
+- #### Clean slate
+
+  Load into a blank simulation, discarding everything. Useful for creating new simulations.
+
+- #### Ghost BCs
+
+  Toggle the availability of `ghost' boundary conditions in the GUI.
+  JSON key: `showGhostBCs`, boolean.
+
+- #### Run on load
+
+  Toggle whether or not the simulation will be running when loaded.
+  JSON key: `runningOnLoad`, boolean.
+
+- #### Blend image
+
+  Toggle whether or not an image will be blended over the simulation.
+  JSON key: `blendImage`, boolean.
+
+- #### Blend amount
+
+  The degree to which the blend image should be blended into the simulation view. 0 corresponds to no blending.
+  JSON key: `blendImageAmount`, float.
+
+- #### Dev
+
+  Tools intended for the development and benchmarking of VisualPDE.
+
+    ***Copy code*** will copy a verbose description of your simulation in JSON form, which is especially useful if you're extending VisualPDE with your own examples. It will base the example on the selected 'parent' preset, which can be useful if you're making multiple slightly different versions of a simulation.
+
+    ***Copy debug*** will copy a selection of configuration information to your clipboard (handy when reporting bugs).
+
+    ***Show stats*** will toggle the display of performance statistics in the lower left corner of the display. We use this to benchmark performance impacts of new features and optimisations.
+    JSON key: `showStats`, boolean.
+
+    ***Antialias*** will toggle the use of antialiasing when displaying the simulation. Antialiasing in VisualPDE smooths out the jagged edges of displayed vector fields, though always causes a slight device-dependent performance loss. This is off by default for mobile devices, but can be overridden. Toggle requires a page reload.
+
+    ***Mixed strings*** will toggle the visibility of verbose strings associated with mixed boundary conditions. These are off by default, but can be useful for advanced cases such as simultaneously imposing Dirichlet and Neumann conditions.
+
+    ***Long url*** will copy a long URL that is guaranteed to not be minified.
+
+    ***Camera*** will enable the user's camera and replace $I_S$ with the live video stream.
+
+    ***Cam delay*** sets the delay (in milliseconds) between frames grabbed from the user's camera.
+
+---
+
+## Writing valid expressions
+
+### Standard syntax
+
+VisualPDE aims to support standard mathematical syntax for addition `+`, subtraction `-`, multiplication `*`, and division `/`, along with the caret notation `^` for exponentiation. Parentheses `()` are also supported for bracketing terms together, and must be used when calling any special functions.
+
+In general, whitespace around binary operators will be ignored, so that `2 * 2` is valid syntax for multiplication, for example. Quantities with subscripts, such as $L_x$ and $L_y$, are written with an underscore, e.g. `L_x` and `L_y`. Boolean expressions (e.g. `x < 0`) can be used as indicator functions by surrounding them with the function `ind` (e.g. `ind(x < 0)`).
+
+Examples of valid syntax include
+
+```
+sin(x) * cos(3*y)
+exp( -(x-L_x)^2 / 10)
+1 + (x^2 + x + 1) / (y^2 + 2*y + 1)
+sin(cosh(tan(2*x+1)))
+```
+
+### Advanced syntax
+
+Some terms in VisualPDE have additional functionality when written with special syntax.
+
+#### Images
+
+By default, images are accessed using `I_T` and `I_S`, with individual channels available by appending R,G,B or A. When a channel is specified, you can access images using coordinates using the syntax `I_TR(x,y)`. Examples include
+
+```
+I_TR(2*x, y)
+I_SB(x-t, y-t)
+I_TG(x, y + sin(u))
+```
+
+#### First derivatives
+
+First derivatives in space, accessed with `u_x`, `u_y`, ..., are computed using a central finite difference discretisation by default. By appending `f` or `b` to the subscript, such as `u_xf`, you can tell VisualPDE to use a forward or a backward difference, respectively. Forward differences sample the solution at increased $x$ (or $y$), whilst backward differences sample at decreased $x$ (or $y$). These specialised schemes can be used in [upwind schemes](https://en.wikipedia.org/wiki/Upwind_scheme) and often reduce numerical artefacts, but at the expense of typically larger numerical error.
+
+Forward and backward differences can also be computed with second-order numerical schemes by appending `2` to the subscript, though in general this will only respect Periodic boundary conditions in the direction of the derivative. This syntax can only be used in the **Edit** section of **Equations**.
+
+### Special functions
+
+Throughout VisualPDE, you can make use of the special functions `sin`, `cos`, `tan`, `exp`, `log`, `sqrt`, `sinh`, `cosh`, `tanh` and `H`, where the latter is a [Heaviside function](https://en.wikipedia.org/wiki/Heaviside_step_function) smoothed over the interval $[-1,1]$ (see the [GLSL reference](https://registry.khronos.org/OpenGL-Refpages/gl4/html/smoothstep.xhtml) for details). All function arguments should be surrounded by parentheses, e.g. `sin(x)`. You can also use `min` and `max` as functions with two arguments, which return the minimum or maximum of their arguments, e.g. `min(u,1)` returns the minimum of $u$ and 1. If you wish to raise the output of a function to a power, you must enclose the function in parentheses, e.g. write `(cos(x))^2`, not `cos(x)^2`. You can also use `mod(a,b)` to compute the remainder of $a$ upon division by $b$ (see the [GLSL reference](https://registry.khronos.org/OpenGL-Refpages/gl4/html/mod.xhtml) for details).
+
+A [bump function](https://en.m.wikipedia.org/wiki/Bump_function) with compact support can be used via the syntax `Bump(X, Y, radius)` (or `Bump(X, radius)` in 1D, with `Y` implicitly set to `L_y/2`), which localises a bump of unit maximum of the given radius at the point $(X, Y)$.
+
+A [bivariate Gaussian function](https://en.wikipedia.org/wiki/Multivariate_normal_distribution) can be used with similar syntax to the bump function: `Gauss(X, Y, s)` is a correctly normalised Gaussian function centred at $(X,Y)$ with standard deviation `s`. The extended syntax `Gauss(X, Y, s_x, s_y)` produces a potentially asymmetric Gaussian with standard deviations `s_x` and `s_y` in the $x$ and $y$ directions, with zero correlation. Correlation can be specified via `Gauss(X, Y, s_x, s_y, r)`, where `r` is the correlation between the two directions.
+
+### Non-local evaluation
+
+Sometimes, you might want direct access to the values of a species at points other than `(x,y)`, such as if you wanted to implement a higher order derivative or a lattice dynamical system. You can access a species (e.g. `u`) at an arbitrary position via the syntax `u[x,y]`, where `x` and `y` can be any valid expression that points to somewhere in the domain. For instance, writing `u[x+dx,y+dy] - u[x-dx,y-dy]` will lookup nearby values of `u`, offset by the step sizes `dx` and `dy`.
+
+The behaviour of non-local evaluations using these expressions is undefined at boundaries (so that Neumann boundary conditions may not behave as expected, for instance), unless periodic boundary conditions are used.
+
+## Interaction via query strings
+
+The simulation can be modified directly using query strings.
+
+For example `https://visualpde.com/sim/?boundaryConditions_1=dirichlet` replaces the default simulation preset first variable's boundary condition to Dirichlet.
+
+This can be chained using standard query string notation using `&`. For instance, `https://visualpde.com/sim/?boundaryConditions_1=dirichlet&timesteppingScheme=RK4` replaces the default simulation preset first variable's boundary condition to Dirichlet *and* the timestepping scheme to Runge-Kutta 4.
+
+### Other parameters
+
+The following parameters are only accessible via query strings.
+
+- #### Simulation Title
+
+  Set the title of the simulation shown at the top left of the interface.
+  JSON key: `simTitle`, string. Default: `Grey-Scott`.
+
+### Known quirks
+
+There are special characters that are reserved in query strings that often overlap with mathematical expressions, these are,
+
+- To denote space we use `+` or `%20`, commonly used to rename species, e.g. `speciesNames=u+v+w`.
+- Conversely, to denote a literal plus symbol `+`, we use `%2B`, common in expressions, e.g. `reactionStr_3=v*w+0.5*w`.
+- Similarly, `%2F` is used for `/`, commonly used as the division operator, e.g. `reactionStr_1=v*w/u`.
+- Lastly, `%26` is used for `&`, but this is not used in VisualPDE.
+
+### Example
+
+The following gives a 3-species system:
+
+`https://visualpde.com/sim/?numSpecies=3&speciesNames=u+v+w&crossDiffusion=true&reactionStr_1=u*(1-u)-u*v&reactionStr_2=u*v-0.5*v-v*w&reactionStr_3=v*w-0.5*w&diffusionStr_1_1=0.1&diffusionStr_2_2=0.1&diffusionStr_3_3=0.1&diffusionStr_1_2=0.5&initCond_1=RAND&initCond_2=0.1&initCond_3=0.01&boundaryConditions_1=dirichlet&dirichletStr_1=0&boundaryConditions_2=neumann&neumannStr_2=0&boundaryConditions_3=periodic&dt=0.001&simTitle=3-species%20example`
+
+Explaination,
+
+- `numSpecies=3`: sets `# species` to 3.
+- `speciesNames=u+v+w`: sets species names to `u`, `v` and `w`. This is essential as these will be used in the expressions that follow.
+- `crossDiffsuin=true`: toggles cross diffusion.
+- `reactionStr_1=u*(1-u)-u*v`: sets `f_u = u*(1-u)-u*v`.
+- `reactionStr_2=u*v-0.5*v-v*w`: sets `f_v = u*v-0.5*v-v*w`.
+- `reactionStr_3=v*w-0.5*w`: sets `f_w = v*w-0.5*w`.
+- `diffusionStr_1_1=0.1`: sets `D_{uu} = 0.1`.
+- `diffusionStr_2_2=0.1`: sets `D_{vv} = 0.1`.
+- `diffusionStr_3_3=0.1`: sets `D_{ww} = 0.1`.
+- `diffusionStr_1_2=0.5`: sets `D_{uv} = 0.5`.
+- `initCond_1=RAND`: sets `u|_{t=0}=RAND`.
+- `initCond_2=0.1`: sets `v|_{t=0}=0.1`.
+- `initCond_3=0.01`: sets `u|_{t=0}=0.01`.
+- `boundaryConditions_1=dirichlet`: sets `u` boundary condition to Dirichlet.
+- `dirichletStr_1=0`: sets `u|_{\partial\Omega}=0`.
+- `boundaryConditions_2=neumann`: sets `v` boundary condition to Neumann.
+- `neumannStr_2=0`: sets `\frac{\partial v}{\partial n}|_{\partial\Omega}=0`.
+- `boundaryConditions_2=periodic`: sets `v` boundary condition to Periodic.
+- `dt=0.001`: sets the timestep size to `0.001`. The default timestep size of `0.1` is too large for this kind of simulation.
+- `simTitle=3-species%20example`: sets the simulation title to `3-species example`.
+
+All other options remain to default values as they are unmodified.
diff --git a/sim/llm-context/misc.md b/sim/llm-context/misc.md
new file mode 100644
index 00000000..f2e8cf0e
--- /dev/null
+++ b/sim/llm-context/misc.md
@@ -0,0 +1,17 @@
+# Misc
+
+## Simulator defaults
+
+### Starting a new simulation
+
+Starting a new simulation will start a Grey-Scott simulation with two species, $u$ and $v$.
+
+### Changing the default variable names
+
+When changing "# species" to a higher number previously set, new species will be named by default to `SPECIES2`, `SPECIES3` and `SPECIES4`, for 2, 3, 4 species, respectively.
+
+You are advised to change them to $v$, $w$, $q$ for simplicity, as other parameters are also defined according to them, e.g. $D_{SPECIES2} = D_v$ and so on.
+
+### Default timestepping
+
+By default, VisualPDE uses the Forward Euler method with $\Delta t = 0.1$, this is likely to be too big for many problems. We recommend reducing the timestep and/or changing the method in this case.
diff --git a/sim/llm-context/quick-start.md b/sim/llm-context/quick-start.md
new file mode 100644
index 00000000..27a831e1
--- /dev/null
+++ b/sim/llm-context/quick-start.md
@@ -0,0 +1,75 @@
+# Quick start
+
+VisualPDE is a web-based set of tools for solving partial differential equations (PDEs) via an interactive, easy-to-use simulation. To get started, try playing with some of the [linear examples](/explore#linear), or read on for some quick tips for using the solver.
+
+### Interacting with the simulation
+Clicking/pressing on the simulation draws values right onto the domain. You can customise exactly what this does under <span class='click_sequence'>Settings (🔧) → **Brush** </span> For example, the default settings in the [heat equation example](/basic-pdes/heat-equation) allow you to paint 'heat' of value 1 onto the domain, which acts like an initial condition for the rest of the simulation. If you have another mouse button to press, you can use it to remove 'heat'.
+
+### The equations panel
+Pressing Equations ($f(x)$) opens up the **equations panel**.
+
+Here you can:
+* See the [equation being simulated](#equations), here $\frac{du}{dt} = \nabla\cdot(D_u\nabla u) + f_u$.
+* Set the named functions in the equations, here $D_u$ and $f_u$, under **Definitions**. These can be functions of any of the unknowns, space, and time (here $u$, $x$, $y$, and $t$), and of any parameters that will be defined further down the panel.
+* Set the value of any extra parameters.
+* Set the [boundary conditions](#boundary-conditions).
+* Set the [initial conditions](#initial-conditions).
+* Set the [number and type of equations](#changing-the-equations-) to be solved.
+
+### Domain shape
+
+The default **domain** for solving PDEs is a 2D rectangle, $\Omega = [0,L_x]\times[0,L_y]$, which fits the size of your browser window or phone screen. Throughout VisualPDE, we use coordinates $x\in[0,L_x]$ and $y\in[0,L_y]$.
+
+You can force the domain to be a square, $\Omega = [0,L]\times[0,L]$, by toggling off <span class='click_sequence'>Settings (🔧) → **Domain** → **Fill screen**</span>
+
+### Boundary conditions
+
+The following **boundary conditions** are available to allow you to set the value of the function, or the value of its derivative, along the boundary $\partial\Omega$ of the domain $\Omega$:
+
+* Periodic
+* [Dirichlet](https://en.wikipedia.org/wiki/Dirichlet_boundary_condition) (e.g. $u|_{\partial\Omega} = 0$)
+* [Neumann](https://en.wikipedia.org/wiki/Neumann_boundary_condition) (e.g. $\frac{du}{dn}|_{\partial\Omega} = 0$)
+* [Robin](https://en.wikipedia.org/wiki/Robin_boundary_condition) (e.g. $(u + \frac{du}{dn})|_{\partial\Omega} = 0$)
+
+You can swap between boundary conditions by choosing <span class='click_sequence'>Equations ($f(x)$) → **Boundary conditions**</span> and selecting from the list for each variable.
+
+### Initial conditions
+You can specify the values to which the unknowns ($u$, $v$, $w$) are initialised when resetting the simulation. These expressions can be functions of $x$, $y$, the special string 'RAND' that assigns a random number in [0,1] to each point in the domain, along with any user-defined parameters and the images $I_S$ and $I_T$ (see the [advanced documentation](/user-guide/advanced-options) for more details). You can also use $L$, $L_x$ and $L_y$.
+
+### Changing the equations
+
+The simplest system VisualPDE can solve is a single PDE,
+
+$$\frac{du}{dt} = \nabla \cdot (D_u \nabla u) + f_u,$$
+
+where $D_u$ and $f_u$ are functions of $u$, $x$, $y$, and $t$ that you can specify.
+
+The most complicated type is a coupled system of PDEs in four unknowns, $u$, $v$, $w$ and $q$:
+
+$$\begin{aligned}
+t_u\frac{du}{dt} &= \nabla \cdot(D_{uu}\nabla u+D_{uv}\nabla v+D_{uw}\nabla w+D_{uq}\nabla q) + f_u,\\
+\text{one of}\left\{\begin{matrix}\displaystyle t_v\frac{dv}{dt} \\ v\end{matrix}\right. & 
+\begin{aligned}
+    &= \nabla \cdot(D_{vu}\nabla u+D_{vv}\nabla v+D_{vw}\nabla w+D_{vq}\nabla q) + f_v \vphantom{\displaystyle t_v\frac{dv}{dt}}, \\
+    &= \nabla \cdot(D_{vu}\nabla u+D_{vw}\nabla w+D_{vq}\nabla q) + f_v,
+\end{aligned}\\
+\text{one of}\left\{\begin{matrix}\displaystyle t_w\frac{dw}{dt} \\ w\end{matrix}\right. & 
+\begin{aligned}
+    &= \nabla \cdot(D_{wu}\nabla u+D_{wv}\nabla v+D_{ww}\nabla w+D_{wq}\nabla q) + f_w \vphantom{\displaystyle t_w\frac{dw}{dt}}, \\
+    &= \nabla \cdot(D_{wu}\nabla u+D_{wv}\nabla v+D_{wq}\nabla q) + f_w,
+\end{aligned}\\
+\text{one of}\left\{\begin{matrix}\displaystyle t_q\frac{dq}{dt} \\ q\end{matrix}\right. & 
+\begin{aligned}
+    &= \nabla \cdot(D_{qu}\nabla u+D_{qv}\nabla v+D_{qw}\nabla w+D_{qq}\nabla q) + f_q \vphantom{\displaystyle t_q\frac{dq}{dt}}, \\
+    &= \nabla \cdot(D_{qu}\nabla u+D_{qv}\nabla v+D_{qw}\nabla w) + f_q,
+\end{aligned}
+\end{aligned}$$
+
+where $D_{uu}, \dots,  D_{qq}$, $f_u, \dots, f_q$ and $t_u, \dots, t_q$ are functions of $u$, $v$, $w$, $q$, $x$, $y$ and $t$ that you can specify.
+
+* You can change the number of unknowns by choosing <span class='click_sequence'>Equations ($f(x)$) → **Advanced options** → **\# species**</span>
+* In systems of multiple unknowns, you can include terms representing cross-diffusion (e.g. $D_{uv}$, $D_{vu}$) by toggling <span class='click_sequence'>Equations ($f(x)$) → **Advanced options** → **Cross diffusion**</span>
+* In systems of multiple unknowns, you can choose between a differential or algebraic equation for some of the species (e.g. '$\partial w/\partial t=$' or '$w=$') by toggling <span class='click_sequence'>Equations ($f(x)$) → **Advanced options** → **Algebraic w** (or **v** or **q**)</span>
+
+### More VisualPDE
+For a comprehensive list of all the options that you can set in VisualPDE, check out the [Advanced documentation](/user-guide/advanced-options), or discover what VisualPDE can solve in our brief [summary](/user-guide/what-can-visualpde-solve).
diff --git a/sim/llm-context/solver.md b/sim/llm-context/solver.md
new file mode 100644
index 00000000..250311ea
--- /dev/null
+++ b/sim/llm-context/solver.md
@@ -0,0 +1,91 @@
+# The VisualPDE solver
+
+VisualPDE aims to be a plug-and-play, browser-based solver and visualiser for a range of PDE systems, with as few barriers to entry as possible. Underneath what we hope is a polished exterior, we've made a number of somewhat technical choices to enable interactive and reliable solution across as many devices as we can support. To explain these choices, and because we like knowing what's going on inside any software we use, we've compiled a short summary of how VisualPDE works under the hood. This is accompanied by a more detailed account in an [open access article](https://doi.org/10.1007/s11538-023-01218-4).
+
+We are always looking for ways to improve and extend VisualPDE, especially ways of reaching a broader audience and new communities. If you have any questions or suggestions about anything related to VisualPDE, we'd love to hear from you at [hello@visualpde.com](mailto:hello@visualpde.com)!
+
+### The equations
+
+VisualPDE can solve a variety of PDE systems posed in 1D or 2D space, many of which are straightforward extensions of the two-species reaction–diffusion system,
+
+$$
+\begin{aligned}
+    \frac{du}{dt} &= \nabla \cdot(D_u\nabla u) + f_u,\\
+    \frac{dv}{dt} &= \nabla \cdot(D_v\nabla v) + f_v.
+\end{aligned}
+$$
+
+Here, $u$ and $v$ are the scalar unknowns that we solve for, $t$ is time, and the divergence and gradient operators are the usual spatial operators in a 2D Euclidean domain. In general, the interaction/kinetic terms ($f_u$ and $f_v$) and the diffusion coefficients ($D_u$ and $D_v$) can each be functions of time, space, and the unknowns, though we often don't explicitly write these potential dependencies. See this [help page](what-can-visualpde-solve) for a comprehensive summary of the types of PDEs that VisualPDE can solve.
+
+### The domain
+
+A PDE problem is not well-posed without specifying a domain (or [boundary conditions](#boundary-conditions)). We typically determine the domain $\Omega$ from the size of your device, fixing the largest side to be of length $L$. You can demand a square domain with aspect ratio 1:1 by toggling <span class='click_sequence'>Settings (🔧) → **Domain** → **Fill screen**</span>
+
+However, not all domains are rectangular. To accommodate this, VisualPDE allows you to specify an arbitrary domain via a [level set](https://en.wikipedia.org/wiki/Level-set_method) or [indicator function](https://en.wikipedia.org/wiki/Indicator_function) approach, under <span class='click_sequence'>Settings (🔧) → **Domain** → **Implicit**</span> To make this work in practice, VisualPDE uses this user-specified function to determine which parts of $\Omega$ to include in the simulation, though this prevents you from specifying boundary conditions that include derivatives (computing normals to user-specified curves is not something we wanted to think about). You can specify a boolean (e.g. $x<0.5$) or a simple expression (e.g. $x-0.5$), where (strict) positivity identifies the interior of the domain. You can even use images in this expression, allowing you to define complicated domains with ease.
+
+### Spatial discretisation
+
+We take $\Omega$ and divide it into a square grid with a user-configurable spacing. Be warned: the only limit to this spacing is your imagination/hardware, so taking a spatial step of $\Delta x=10^{-9}$ is unlikely to be a good idea... We use this grid in a [finite difference](https://en.wikipedia.org/wiki/Finite_difference_method) scheme to approximate the spatial operators in our system.
+
+If $D_u$ and $D_v$ were simply constants, as they were when we began the development of VisualPDE, a standard [central differences](https://en.wikipedia.org/wiki/Discrete_Laplace_operator) discretisation of the resulting Laplacian $\nabla^2$ would suffice. However, as these coefficients generally vary in space, we employ a similar but necessarily more complex scheme. Explicitly, using $\nabla \cdot(D_u\nabla u)$ as an example and limiting ourselves to 1D for brevity, we approximate
+
+$$\textstyle \nabla \cdot(D_u\nabla u) \approx \frac{D_u(x-\Delta x)[u(x-\Delta x) - u(x)] + D_u(x)[u(x-\Delta x) - 2u(x) + u(x+\Delta x)] + D_u(x+\Delta x)[u(x+\Delta x) - u(x)]}{2\,\Delta x^2}$$
+
+at a point $x$, where we've omitted any dependence of any quantities on anything other than space. Notably, this is just a standard central differences scheme if $D_u$ is constant. Adding the above to itself with $x$ replaced with $y$ gives the 2D discretisation, and adding in additional terms is simple by linearity.
+
+### Timestepping
+
+With space discretised as above, we are faced with a large system of coupled ordinary differential equations to solve, which represent the evolution of the unknowns at each discrete gridpoint of the spatial domain. VisualPDE implements four popular schemes for timestepping: [Forward Euler](https://en.wikipedia.org/wiki/Euler_method), two-step [Adams-Bashforth](https://en.wikipedia.org/wiki/Linear_multistep_method#Two-step_Adams–Bashforth), the [Midpoint Method](https://en.wikipedia.org/wiki/Midpoint_method) and the four-step [Runge-Kutta](https://en.wikipedia.org/wiki/Runge–Kutta_methods) method (commonly known as 'RK4'). These solvers each have their strengths, with Forward Euler being the least computationally demanding while RK4 offers superior accuracy and stability at the cost of doing more calculations each timestep. The following Forward Euler scheme is the default in many of the examples on the site:
+
+$$ \frac{du}{dt} \approx \frac{u(t+\Delta t) - u(t)}{\Delta t}$$
+
+for timestep $\Delta t$. This approach is far from state-of-the-art, but it is straightforward and intuitive to implement on massively parallel computing hardware (more on that [later](#browser)).
+
+In practice, its simplicity can lead to some problems, with [numerical instability](https://en.wikipedia.org/wiki/Euler_method#Numerical_stability) being perhaps the most pathological. Loosely speaking, the interaction between the forward Euler scheme and our spatial discretisation can lead to numerical artefacts ruining the solution, which typically occur when the ratio $D \, \Delta t / \Delta x^2$ is too small, where $D$ is any of the diffusion coefficients in the problem. Both the Midpoint Method and RK4 improve upon the stability of Forward Euler, whilst the Adams-Bashforth scheme is generally less stable but more accurate than Forward Euler.
+
+VisualPDE will try to tell you when it's fallen foul of stability issues (we periodically check for 'NaN' or $\pm\infty$ in the solution), at which point you might want to try reducing $\Delta t$, trying out a different solver, reducing the diffusion coefficients in your problem, or increasing $\Delta x$ (we recommend trying each of these things in this order).
+Despite each of our solvers having their limitations, these schemes have enabled VisualPDE to efficiently solve every system that we've thrown at it, though some tuning of the timestep can be necessary in extreme cases. If you have any tips for implementing alternative schemes (especially anything implicit), we'd love to hear from you!
+
+### Boundary conditions
+
+VisualPDE implements four types of boundary condition: periodic, [Dirichlet](https://en.wikipedia.org/wiki/Dirichlet_boundary_condition), [Neumann](https://en.wikipedia.org/wiki/Neumann_boundary_condition), and [Robin](https://en.wikipedia.org/wiki/Robin_boundary_condition). As each of these are slightly different in character, we briefly describe the general form of each condition that can be used in VisualPDE, along with notes on how this is enforced in the simulation. We pose these conditions as if they correspond to the scalar [heat equation](/\_basic-pdes/heat-equation) in 2D.
+
+#### Periodic
+
+Periodic boundary conditions, where $u(x+L_x,y)=u(x,y)$ etc, are the simplest to implement. If our spatial discretisation is attempting to sample a function outside of $\Omega$, we simply map its argument back to $\Omega$ by shifting by $L_x$ or $L_y$.
+
+#### Dirichlet
+
+Dirichlet boundary conditions take the form $u|_{\partial\Omega} = a(x,y,t)$ for a user-specified function $a$. These are easy to enforce, as we simply override our usual timestepping for any nodes on the boundary of the domain. If the user has chosen to use an implicitly defined domain, we assign $a(x,y,t)$ to all points outside of $\Omega$.
+
+#### Neumann
+
+Neumann boundary conditions are specified as $\frac{du}{dn}|_{\partial\Omega} = a(x,y,t)$ for a user-specified function $a$, where $\frac{du}{dn}$ denotes a derivative in the direction of the (outward-pointing) normal to the boundary. Implementing a Neumann boundary condition is done via so-called **ghost nodes** in our discretisation. For instance, enforcing $\frac{du}{dn}|_{\partial\Omega} = 0$ at the left-most $x$ boundary of a rectangular domain is achieved in practice by taking
+
+$$\textstyle u(x-\Delta x,y) = u(x+\Delta x,y)$$
+
+in the [finite difference operator](#spatial-discretisation) described above.
+
+#### Robin
+
+Robin boundary conditions are a natural combination of Dirichlet and Neumann conditions, which we pose in the form of a generalised Neumann condition $\frac{du}{dn}|_{\partial\Omega} = a(u,x,y,t)$, where the right-hand side can now depend on $u$ (and any other unknown in multi-species systems). These conditions are also implemented with ghost nodes. For example, enforcing $\frac{du}{dn}|_{\partial\Omega} = u|_{\partial\Omega}$ at the leftmost $x$ boundary of a rectangular domain is achieved in practice by taking
+
+$$\textstyle u(x-\Delta x,y) = u(x+\Delta x,y) + 2 u(x,y)\,\Delta x$$
+
+in the [finite difference operator](#spatial-discretisation) described above, approximating the derivative at the boundary with a simple central difference.
+
+#### Combined
+
+VisualPDE also allows you to specify different boundary conditions on different parts of the boundary. Doing this requires some special syntax, as detailed in the [advanced documentation](/user-guide/advanced-options#boundary-conditions).
+
+### Doing this in your browser, quickly
+
+Solving PDEs is hard. To solve them in real time in your browser, VisualPDE gives all the hard work to the graphics chip (GPU) on your device, making use of [WebGL](https://en.wikipedia.org/wiki/WebGL) and a low-level shader language called [GLSL](https://en.wikipedia.org/wiki/OpenGL_Shading_Language).
+
+Every time your browser requests a frame from VisualPDE (which might be up to 60 times per second), some [JavaScript](https://en.wikipedia.org/wiki/JavaScript) organises the solving of the discretised equations, displaying the solution, and incorporating anything you've drawn, which all happen on the GPU. Each frame, we typically perform hundreds of timesteps to give you a smooth experience, mitigating many of the limitations of our [timestepping schemes](#timestepping). If you're interested in the finest details of the implementation, the source code for the entire site is freely available to view, reuse, and repurpose on [GitHub](https://github.com/Pecnut/visual-pde).
+
+### Accuracy and precision
+
+VisualPDE hopes to be as accurate as possible whilst providing a responsive, visual, portable platform for solving PDEs. In most systems, the timestep, timestepping scheme and spatial discretisation will be the main source of any errors, as you'd expect from finite-difference discretisations of PDEs. Naturally, smaller timesteps, higher order timestepping schemes and refined spatial discretisations will often improve the accuracy of the solution, but each will incur additional computational costs. With VisualPDE, you can choose the balance that works best for you (and potentially your audience).
+
+A more subtle limit on the accuracy of VisualPDE is our use of single-precision arithmetic, something we've inherited from the libraries we use and the capabilities of modern hardware. We'll be moving to double-precision computation as soon as we can to maximise the accuracy of VisualPDE. In the meantime, rest assured that single-precision arithmetic is sufficient for capturing a wide range of phenomena, including those explored in our examples, but keep this caveat in mind when seeking to find precise quantitative answers to PDE problems.
diff --git a/sim/llm-context/what-can-visualpde-solve.md b/sim/llm-context/what-can-visualpde-solve.md
new file mode 100644
index 00000000..851499ba
--- /dev/null
+++ b/sim/llm-context/what-can-visualpde-solve.md
@@ -0,0 +1,44 @@
+# What can VisualPDE solve?
+
+VisualPDE solves systems of PDEs that look like generalised reaction--diffusion equations. It can do this in 1D or 2D.
+
+The simplest type of system is just a single PDE in a single unknown, $u$,
+
+$$\frac{du}{dt} = \nabla \cdot (D_u \nabla u) + f_u,$$
+
+where $D_u$ and $f_u$ are functions of $u$, $t$, and space that you can choose. For example, if $f_u=0$ and $D_u$ is a constant, you have [the heat equation](/basic-pdes/heat-equation). 
+
+The most complicated type is a coupled system of PDEs in four unknowns, $u$, $v$, $w$ and $q$:
+
+$$\begin{aligned}
+t_u\frac{du}{dt} &= \nabla \cdot(D_{uu}\nabla u+D_{uv}\nabla v+D_{uw}\nabla w+D_{uq}\nabla q) + f_u,\\
+\text{one of}\left\{\begin{matrix}\displaystyle t_v\frac{dv}{dt} \\ v\end{matrix}\right. & 
+\begin{aligned}
+    &= \nabla \cdot(D_{vu}\nabla u+D_{vv}\nabla v+D_{vw}\nabla w+D_{vq}\nabla q) + f_v \vphantom{\displaystyle t_v\frac{dv}{dt}}, \\
+    &= \nabla \cdot(D_{vu}\nabla u+D_{vw}\nabla w+D_{vq}\nabla q) + f_v,
+\end{aligned}\\
+\text{one of}\left\{\begin{matrix}\displaystyle t_w\frac{dw}{dt} \\ w\end{matrix}\right. & 
+\begin{aligned}
+    &= \nabla \cdot(D_{wu}\nabla u+D_{wv}\nabla v+D_{ww}\nabla w+D_{wq}\nabla q) + f_w \vphantom{\displaystyle t_w\frac{dw}{dt}}, \\
+    &= \nabla \cdot(D_{wu}\nabla u+D_{wv}\nabla v+D_{wq}\nabla q) + f_w,
+\end{aligned}\\
+\text{one of}\left\{\begin{matrix}\displaystyle t_q\frac{dq}{dt} \\ q\end{matrix}\right. & 
+\begin{aligned}
+    &= \nabla \cdot(D_{qu}\nabla u+D_{qv}\nabla v+D_{qw}\nabla w+D_{qq}\nabla q) + f_q \vphantom{\displaystyle t_q\frac{dq}{dt}}, \\
+    &= \nabla \cdot(D_{qu}\nabla u+D_{qv}\nabla v+D_{qw}\nabla w) + f_q,
+\end{aligned}
+\end{aligned}$$
+
+where the diffusion coefficients ($D_{uu}$ etc.), the timescales ($t_u$ etc.) and the interaction/kinetic terms ($f$, $g$, $h$, $j$) can depend on the unknowns, space, and time. In matrix form, we can summarise this by saying we solve systems of the form
+
+$$M \frac{du}{dt} = \nabla\cdot(D\nabla u) + f,$$
+
+where
+
+* $u$ is a vector of one, two, three or four unknowns,
+* $M$ is a diagonal matrix with potentially some zeros on the diagonal; you might know this as a 'mass matrix',
+* $D$ is a possibly non-constant matrix that may contain zeros; you might know this as a 'diffusion tensor',
+* $f$ is a vector of one, two, three or four components that contains our interaction or kinetic terms.
+
+VisualPDE allows you to easily change the [number of components](quick-start#equations-panel) and the [boundary conditions](quick-start#boundary-conditions). You can set initial conditions just by clicking the screen.
+
diff --git a/sim/scripts/RD/main.js b/sim/scripts/RD/main.js
index dcd6186a..b054569e 100644
--- a/sim/scripts/RD/main.js
+++ b/sim/scripts/RD/main.js
@@ -91,6 +91,7 @@ import {
 import { closestMatch } from "../../../assets/js/closest-match.js";
 import { Stats } from "../stats.min.js";
 import { createWelcomeTour } from "./tours.js";
+import { GoogleGenAI } from "@google/genai";
 
 if (expandingOptionsInProgress) {
   let checkIfOptionsLoaded = setInterval(() => {
@@ -962,6 +963,72 @@ async function VisualPDE(url) {
     );
   }
 
+  const aiBtn = document.getElementById("ai-generate-btn");
+  const aiInput = document.getElementById("ai-prompt-input");
+  const aiStatus = document.getElementById("ai-status");
+  const ai = new GoogleGenAI({ apiKey: "" });
+  const manifest = {
+      "what-can-visualpde-solve.md": { title: "What can VisualPDE solve?", url: "https://visualpde.com/user-guide/what-can-visualpde-solve.html" },
+      "FAQ.md": { title: "FAQs", url: "https://visualpde.com/user-guide/FAQ.html" },
+      "quick-start.md": { title: "Quick Start", url: "https://visualpde.com/user-guide/quick-start.html" },
+      "advanced-options.md": { title: "Complete Documentation", url: "https://visualpde.com/user-guide/advanced-options.html" },
+      "solver.md": { title: "The VisualPDE solver", url: "https://visualpde.com/user-guide/solver.html" },
+      "misc.md": { title: "Misc information", url: "N/A" }
+  };
+
+  async function loadDocs() {
+      let combinedDocs = "OFFICIAL DOCUMENTATION CONTENT:\n";
+      
+      for (const [filename, info] of Object.entries(manifest)) {
+          try {
+              const response = await fetch(`/sim/llm-context/${filename}`);
+              const text = await response.text();
+              
+              combinedDocs += `\n\n--- SECTION: ${info.title} ---\n`;
+              combinedDocs += `SOURCE URL: ${info.url}\n`;
+              combinedDocs += text;
+          } catch (err) {
+              console.error(`Could not load ${filename}:`, err);
+          }
+      }
+      return combinedDocs;
+  }
+
+  const SYSTEM_INSTRUCTION = `You are a helpful assistant proficient in reaction-diffusion equations and the official support bot for VisualPDE (visualpde.com).
+  Your job is, according to user query, to ONLY provide simulations raw links by using query strings. 
+  Use the provided documentation as your main reference point. If the query falls outside VisualPDE's capabilities, return 'visualpde.com/sim'.`
+
+  async function generateSim() {
+      const userQuery = document.getElementById("ai-prompt-input").value;
+      const status = document.getElementById("ai-status");
+      
+      status.innerText = "Thinking..";
+
+      const docsContext = await loadDocs();
+
+      const contents = [
+          { text: `DOCUMENTATION CONTEXT:\n${docsContext}` },
+          { text: `USER QUESTION: ${userQuery}` }
+      ];
+
+      try {
+          const result = await ai.models.generateContent({
+            model: "gemini-3-flash-preview",
+            contents: contents,
+            config: {
+              systemInstruction: SYSTEM_INSTRUCTION,
+            },
+          });
+          const link = result.text;
+          status.innerHTML = `Success! <a href="${link}" style="color: cyan;">Link Generated</a>`;
+      } catch (error) {
+          status.innerText = "Error connecting to Gemini.";
+          console.error(error);
+      }
+  }
+
+  aiBtn.addEventListener("click", generateSim);
+
   //---------------
 
   // Initialise all aspects of the site, including both the simulation and the GUI.