diff --git a/lib/StochasticDiffEq/docs/Project.toml b/lib/StochasticDiffEq/docs/Project.toml
new file mode 100644
index 00000000000..c96d1905710
--- /dev/null
+++ b/lib/StochasticDiffEq/docs/Project.toml
@@ -0,0 +1,40 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
+StochasticDiffEqCore = "19c5a474-6cd1-4a5f-be79-46dc34e54d7f"
+StochasticDiffEqHighOrder = "0520c28c-50fd-4d16-9c96-902fc80b3bab"
+StochasticDiffEqIIF = "ebf54054-c36b-4494-9aba-657e36df524f"
+StochasticDiffEqImplicit = "5080b986-4c76-4669-b5dc-373a41579d5b"
+StochasticDiffEqLeaping = "aefaaa88-39f2-4e89-b162-d61b7e8cc81b"
+StochasticDiffEqLowOrder = "d15fe365-ce7f-450a-828a-7985bd5b681b"
+StochasticDiffEqMilstein = "8c95a807-c8e7-4581-8419-890d001af53e"
+StochasticDiffEqROCK = "db241ea8-0e6b-4abc-8f2d-1adff2294fd9"
+StochasticDiffEqRODE = "49714585-0aa1-4f53-b1e6-a9b8c0d5e03f"
+StochasticDiffEqWeak = "af2a2fcd-1c36-4cbe-a6d0-5afda784a085"
+
+[sources]
+StochasticDiffEq = {path = ".."}
+StochasticDiffEqCore = {path = "../../StochasticDiffEqCore"}
+StochasticDiffEqHighOrder = {path = "../../StochasticDiffEqHighOrder"}
+StochasticDiffEqIIF = {path = "../../StochasticDiffEqIIF"}
+StochasticDiffEqImplicit = {path = "../../StochasticDiffEqImplicit"}
+StochasticDiffEqLeaping = {path = "../../StochasticDiffEqLeaping"}
+StochasticDiffEqLowOrder = {path = "../../StochasticDiffEqLowOrder"}
+StochasticDiffEqMilstein = {path = "../../StochasticDiffEqMilstein"}
+StochasticDiffEqROCK = {path = "../../StochasticDiffEqROCK"}
+StochasticDiffEqRODE = {path = "../../StochasticDiffEqRODE"}
+StochasticDiffEqWeak = {path = "../../StochasticDiffEqWeak"}
+
+[compat]
+Documenter = "1"
+StochasticDiffEq = "7"
+StochasticDiffEqCore = "2"
+StochasticDiffEqHighOrder = "2"
+StochasticDiffEqIIF = "2"
+StochasticDiffEqImplicit = "2"
+StochasticDiffEqLeaping = "2"
+StochasticDiffEqLowOrder = "2"
+StochasticDiffEqMilstein = "2"
+StochasticDiffEqROCK = "2"
+StochasticDiffEqRODE = "2"
+StochasticDiffEqWeak = "2"
diff --git a/lib/StochasticDiffEq/docs/make.jl b/lib/StochasticDiffEq/docs/make.jl
new file mode 100644
index 00000000000..41d0d1927d2
--- /dev/null
+++ b/lib/StochasticDiffEq/docs/make.jl
@@ -0,0 +1,35 @@
+using Documenter
+using StochasticDiffEq
+# StochasticDiffEq @reexports its solver subpackages, but the docstrings the @docs
+# blocks reference live in those subpackages, so they must be in `modules` for
+# Documenter to find them.
+using StochasticDiffEqCore, StochasticDiffEqHighOrder, StochasticDiffEqIIF,
+    StochasticDiffEqImplicit, StochasticDiffEqLeaping, StochasticDiffEqLowOrder,
+    StochasticDiffEqMilstein, StochasticDiffEqROCK, StochasticDiffEqRODE,
+    StochasticDiffEqWeak
+
+cp(joinpath(@__DIR__, "Project.toml"), joinpath(@__DIR__, "src", "assets", "Project.toml"),
+    force = true)
+
+# Keep pages.jl separate so DiffEqDocs.jl can include it when aggregating these docs.
+include("pages.jl")
+
+makedocs(
+    sitename = "StochasticDiffEq.jl",
+    authors = "Chris Rackauckas et al.",
+    clean = true,
+    doctest = false,
+    modules = [StochasticDiffEq, StochasticDiffEqCore, StochasticDiffEqHighOrder,
+        StochasticDiffEqIIF, StochasticDiffEqImplicit, StochasticDiffEqLeaping,
+        StochasticDiffEqLowOrder, StochasticDiffEqMilstein, StochasticDiffEqROCK,
+        StochasticDiffEqRODE, StochasticDiffEqWeak],
+    warnonly = [:docs_block, :missing_docs, :eval_block],
+    format = Documenter.HTML(
+        assets = ["assets/favicon.ico"],
+        canonical = "https://docs.sciml.ai/StochasticDiffEq/stable/"
+    ),
+    pages = pages
+)
+
+# Note: these pages are aggregated into the unified SciML docs by DiffEqDocs.jl
+# (which copies docs/src + docs/pages.jl); deployment happens there.
diff --git a/lib/StochasticDiffEq/docs/pages.jl b/lib/StochasticDiffEq/docs/pages.jl
new file mode 100644
index 00000000000..26381b06754
--- /dev/null
+++ b/lib/StochasticDiffEq/docs/pages.jl
@@ -0,0 +1,21 @@
+pages = [
+    "StochasticDiffEq.jl: SDE solvers and utilities" => "index.md",
+    "Usage" => "usage.md",
+    "Nonstiff Solvers" => [
+        "nonstiff/basic_methods.md",
+        "nonstiff/sra_sri_methods.md",
+        "nonstiff/high_weak_order.md",
+        "nonstiff/commutative_noise.md",
+    ],
+    "Stiff Solvers" => [
+        "stiff/implicit_methods.md",
+        "stiff/split_step_methods.md",
+        "stiff/stabilized_methods.md",
+    ],
+    "Jump Diffusion" => [
+        "jumpdiffusion/tau_leaping.md",
+    ],
+    "Misc Solvers" => [
+        "misc.md",
+    ],
+]
diff --git a/lib/StochasticDiffEq/docs/src/assets/Project.toml b/lib/StochasticDiffEq/docs/src/assets/Project.toml
new file mode 100644
index 00000000000..da1e5c14ff3
--- /dev/null
+++ b/lib/StochasticDiffEq/docs/src/assets/Project.toml
@@ -0,0 +1,6 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
+
+[compat]
+Documenter = "0.27, 1"
\ No newline at end of file
diff --git a/lib/StochasticDiffEq/docs/src/assets/favicon.ico b/lib/StochasticDiffEq/docs/src/assets/favicon.ico
new file mode 100644
index 00000000000..3c6bd470373
Binary files /dev/null and b/lib/StochasticDiffEq/docs/src/assets/favicon.ico differ
diff --git a/lib/StochasticDiffEq/docs/src/index.md b/lib/StochasticDiffEq/docs/src/index.md
new file mode 100644
index 00000000000..18cf5c5cb7d
--- /dev/null
+++ b/lib/StochasticDiffEq/docs/src/index.md
@@ -0,0 +1,81 @@
+# StochasticDiffEq.jl
+
+StochasticDiffEq.jl is a component package in the DifferentialEquations ecosystem for solving stochastic differential equations (SDEs) and random ordinary differential equations (RODEs). It provides a comprehensive suite of high-performance numerical methods for stochastic problems.
+
+## Installation
+
+To install StochasticDiffEq.jl, use the Julia package manager:
+
+```julia
+using Pkg
+Pkg.add("StochasticDiffEq")
+```
+
+## Quick Start
+
+Here's a simple example of solving an SDE:
+
+```julia
+using StochasticDiffEq
+
+# Define the drift and diffusion functions
+function drift!(du, u, p, t)
+    du[1] = 1.01 * u[1]
+end
+
+function diffusion!(du, u, p, t)
+    du[1] = 0.87 * u[1]
+end
+
+# Initial condition and time span
+u0 = [0.5]
+tspan = (0.0, 1.0)
+
+# Define the SDE problem
+prob = SDEProblem(drift!, diffusion!, u0, tspan)
+
+# Solve using the default algorithm
+sol = solve(prob)
+```
+
+## Solver Categories
+
+StochasticDiffEq.jl provides several categories of solvers optimized for different types of problems:
+
+### Nonstiff Solvers
+
+  - **Basic Methods**: Euler-Maruyama, Heun methods
+  - **SRA/SRI Methods**: High-order adaptive methods (SOSRI, SOSRA)
+  - **High Weak Order**: Methods optimized for weak convergence (DRI1)
+  - **Commutative Noise**: Specialized methods for commuting noise terms
+
+### Stiff Solvers
+
+  - **Implicit Methods**: Drift-implicit methods for stiff problems
+  - **Split-Step Methods**: Methods handling stiffness in diffusion
+  - **Stabilized Methods**: SROCK-type methods for parabolic PDEs
+
+### Jump-Diffusion
+
+  - **Tau-Leaping**: Methods for jump-diffusion processes
+
+## Recommended Methods
+
+For most users, we recommend starting with these methods:
+
+  - **General Purpose**: `SOSRI()` - Excellent for diagonal/scalar Itô SDEs
+  - **Additive Noise**: `SOSRA()` - Optimal for problems with additive noise
+  - **Stiff Problems**: `SKenCarp()` - Best for stiff problems with additive noise
+  - **Commutative Noise**: `RKMilCommute()` - For multi-dimensional commutative noise
+  - **High Efficiency**: `EM()` - When computational speed is most important
+
+## Advanced Features
+
+  - Adaptive time stepping with sophisticated error control
+  - Support for all noise types (diagonal, non-diagonal, additive, scalar)
+  - Both Itô and Stratonovich interpretations
+  - Integration with the broader DifferentialEquations.jl ecosystem
+  - GPU compatibility for high-performance computing
+  - Extensive callback and event handling capabilities
+
+See the individual solver pages for detailed information about each method's properties, when to use them, and their theoretical foundations.
diff --git a/lib/StochasticDiffEq/docs/src/jumpdiffusion/tau_leaping.md b/lib/StochasticDiffEq/docs/src/jumpdiffusion/tau_leaping.md
new file mode 100644
index 00000000000..3cac0792658
--- /dev/null
+++ b/lib/StochasticDiffEq/docs/src/jumpdiffusion/tau_leaping.md
@@ -0,0 +1,141 @@
+# Tau-Leaping Methods for Jump-Diffusion
+
+Tau-leaping methods approximate jump processes by "leaping" over multiple potential jump events in a single time step. These methods are essential for efficiently simulating systems with jump-diffusion processes.
+
+## Tau-Leaping Methods
+
+### TauLeaping - Basic Tau-Leaping
+
+```@docs
+TauLeaping
+```
+
+### CaoTauLeaping - Cao's Tau-Leaping
+
+```@docs
+CaoTauLeaping
+```
+
+## Understanding Jump-Diffusion Processes
+
+Jump-diffusion processes combine:
+
+ 1. **Continuous diffusion**: Standard Brownian motion terms
+ 2. **Jump processes**: Discontinuous jumps at random times
+
+General form:
+
+```
+dX = μ(X,t)dt + σ(X,t)dW + ∫ h(X-,z)Ñ(dt,dz)
+```
+
+Where:
+
+  - μ(X,t)dt: Drift term
+  - σ(X,t)dW: Diffusion term
+  - Ñ(dt,dz): Compensated random measure (jumps)
+
+## When to Use Tau-Leaping
+
+**Appropriate for:**
+
+  - Systems with many small jumps
+  - When exact jump simulation is computationally prohibitive
+  - Chemical reaction networks
+  - Population models with birth-death processes
+  - Financial models with rare events
+
+**Not appropriate for:**
+
+  - Systems dominated by large, infrequent jumps
+  - When exact jump timing is critical
+  - Small systems where exact methods are feasible
+
+## Method Characteristics
+
+### TauLeaping:
+
+  - Basic tau-leaping approximation
+  - Fixed tau approach
+  - Good for initial exploration
+
+### CaoTauLeaping:
+
+  - Adaptive tau selection
+  - More sophisticated error control
+  - Better for production simulations
+
+## Configuration
+
+Tau-leaping methods require:
+
+ 1. **Jump rate functions**: λ(X,t) for each reaction/jump type
+ 2. **Jump effects**: How state changes with each jump
+ 3. **Tau selection**: Time step size strategy
+
+```julia
+# Basic setup
+prob = JumpProblem(base_problem, aggregator, jumps...)
+sol = solve(prob, TauLeaping())
+
+# With adaptive tau
+sol = solve(prob, CaoTauLeaping(), tau_tol = 0.01)
+```
+
+## Accuracy Considerations
+
+**Tau-leaping approximation quality depends on:**
+
+  - Jump frequency vs. tau size
+  - State change magnitude per jump
+  - System stiffness
+  - Error tolerance requirements
+
+**Rule of thumb:** Tau should be small enough that jump rates don't change significantly over [t, t+tau].
+
+## Alternative Approaches
+
+**If tau-leaping is inadequate:**
+
+ 1. **Exact methods**: Gillespie algorithm for small systems
+ 2. **Hybrid methods**: Combine exact and approximate regions
+ 3. **Moment closure**: For statistical properties only
+ 4. **Piecewise deterministic**: For systems with rare jumps
+
+## Performance Tips
+
+ 1. **Vectorize jump computations** when possible
+ 2. **Use sparse representations** for large systems
+ 3. **Tune tau carefully** - too large gives poor accuracy, too small is inefficient
+ 4. **Monitor jump frequencies** to validate approximation
+
+## Integration with DifferentialEquations.jl
+
+```julia
+using DifferentialEquations, StochasticDiffEq
+
+# Define base SDE
+function drift!(du, u, p, t)
+    # Continuous drift
+end
+
+function diffusion!(du, u, p, t)
+    # Continuous diffusion
+end
+
+# Define jumps
+jump1 = ConstantRateJump(rate1, affect1!)
+jump2 = VariableRateJump(rate2, affect2!)
+
+# Combine into jump-diffusion problem
+sde_prob = SDEProblem(drift!, diffusion!, u0, tspan)
+jump_prob = JumpProblem(sde_prob, Direct(), jump1, jump2)
+
+# Solve with tau-leaping
+sol = solve(jump_prob, TauLeaping())
+```
+
+## References
+
+  - Gillespie, D.T., "Approximate accelerated stochastic simulation of chemically reacting systems"
+  - Cao, Y., Gillespie, D.T., Petzold, L.R., "Efficient step size selection for the tau-leaping method"
diff --git a/lib/StochasticDiffEq/docs/src/misc.md b/lib/StochasticDiffEq/docs/src/misc.md
new file mode 100644
index 00000000000..9ec41444fd7
--- /dev/null
+++ b/lib/StochasticDiffEq/docs/src/misc.md
@@ -0,0 +1,105 @@
+# Miscellaneous Methods
+
+This page covers specialized methods for particular types of problems or applications.
+
+## Composite Algorithms
+
+### StochasticCompositeAlgorithm - Multi-Method Solving
+
+```@docs
+StochasticCompositeAlgorithm
+```
+
+## RODE Methods (Random ODEs)
+
+### RandomEM - Random Euler Method
+
+```@docs
+RandomEM
+```
+
+### RandomHeun - Random Heun Method
+
+```@docs
+RandomHeun
+```
+
+### RandomTamedEM - Tamed Random Euler
+
+```@docs
+RandomTamedEM
+```
+
+## Langevin Dynamics
+
+### BAOAB - Langevin Integrator
+
+```@docs
+BAOAB
+```
+
+## Predictor-Corrector Methods
+
+### PCEuler - Predictor-Corrector Euler
+
+```@docs
+PCEuler
+```
+
+## Integro-Integral-Form (IIF) Methods
+
+### IIF1M, IIF2M, IIF1Mil - IIF Methods
+
+```@docs
+IIF1M
+```
+
+```@docs
+IIF2M
+```
+
+```@docs
+IIF1Mil
+```
+
+## Simplified Methods
+
+### SimplifiedEM - Simplified Euler-Maruyama
+
+```@docs
+SimplifiedEM
+```
+
+## When to Use Miscellaneous Methods
+
+### StochasticCompositeAlgorithm:
+
+  - When problem characteristics change during integration
+  - Combining methods for different regimes
+  - Automatic method switching based on conditions
+
+### RODE Methods:
+
+  - Random ordinary differential equations
+  - Problems with random parameters but no Brownian motion
+  - Uncertainty quantification applications
+
+### BAOAB:
+
+  - Molecular dynamics simulations
+  - Langevin equations with specific structure
+  - When preserving equilibrium distributions is important
+
+### IIF Methods:
+
+  - Semi-linear problems with stiff linear parts
+  - Problems amenable to integrating factor techniques
+  - When exponential integrators are appropriate
+
+### PCEuler:
+
+  - Problems requiring specific drift-diffusion coupling
+  - When analytical ggprime function is available
+  - Specialized predictor-corrector applications
+
+These methods serve specific niches in stochastic computation and may be optimal for particular problem structures.
diff --git a/lib/StochasticDiffEq/docs/src/nonstiff/basic_methods.md b/lib/StochasticDiffEq/docs/src/nonstiff/basic_methods.md
new file mode 100644
index 00000000000..4d2526d23d9
--- /dev/null
+++ b/lib/StochasticDiffEq/docs/src/nonstiff/basic_methods.md
@@ -0,0 +1,79 @@
+# Basic Nonstiff Methods
+
+This page covers the fundamental explicit methods for solving SDEs. These methods are suitable for non-stiff problems and provide the foundation for more advanced algorithms.
+
+## Euler-Maruyama Methods
+
+### EM - Euler-Maruyama
+
+```@docs
+EM
+```
+
+### EulerHeun - Euler-Heun
+
+```@docs
+EulerHeun
+```
+
+### LambaEM - Adaptive Euler-Maruyama
+
+```@docs
+LambaEM
+```
+
+### LambaEulerHeun - Adaptive Euler-Heun
+
+```@docs
+LambaEulerHeun
+```
+
+## Milstein Methods
+
+### RKMil - Runge-Kutta Milstein
+
+```@docs
+RKMil
+```
+
+## Split Methods
+
+### SplitEM - Split Euler-Maruyama
+
+```@docs
+SplitEM
+```
+
+## When to Use Basic Methods
+
+**Use EM when:**
+
+  - Computational efficiency is most important
+  - Problem is not stiff
+  - Any noise type (most flexible)
+  - Simple implementation needed
+
+**Use EulerHeun when:**
+
+  - Working in Stratonovich interpretation
+  - Need slightly better accuracy than EM
+  - Problem has non-commutative noise
+
+**Use LambaEM/LambaEulerHeun when:**
+
+  - Want adaptive time stepping with basic methods
+  - Need error control but not high accuracy
+  - Good balance of simplicity and adaptivity
+
+**Use RKMil when:**
+
+  - Higher accuracy than Euler methods
+  - Problem has diagonal or scalar noise
+  - Strong order 1.0 convergence required
+
+These methods form the foundation of stochastic numerical analysis. While higher-order methods often provide better performance, the basic methods are essential for:
+
+  - Initial testing and prototyping
+  - Problems where simplicity is preferred
+  - Educational purposes
+  - Fallback options when advanced methods fail
diff --git a/lib/StochasticDiffEq/docs/src/nonstiff/commutative_noise.md b/lib/StochasticDiffEq/docs/src/nonstiff/commutative_noise.md
new file mode 100644
index 00000000000..380fdefb357
--- /dev/null
+++ b/lib/StochasticDiffEq/docs/src/nonstiff/commutative_noise.md
@@ -0,0 +1,125 @@
+# Commutative Noise Methods
+
+When multiple noise sources satisfy commutativity conditions, specialized methods can achieve higher accuracy and efficiency compared to general methods. These methods avoid expensive Lévy area computations while maintaining high order.
+
+## Recommended Commutative Noise Methods
+
+### RKMilCommute - Runge-Kutta Milstein for Commutative Noise
+
+```@docs
+RKMilCommute
+```
+
+### RKMilGeneral - General Milstein for Non-commutative Noise
+
+```@docs
+RKMilGeneral
+```
+
+## Three-Stage Milstein Methods
+
+### WangLi3SMil Family - Fixed Step Milstein Methods
+
+```@docs
+WangLi3SMil_A
+```
+
+```@docs
+WangLi3SMil_B
+```
+
+```@docs
+WangLi3SMil_C
+```
+
+```@docs
+WangLi3SMil_D
+```
+
+```@docs
+WangLi3SMil_E
+```
+
+```@docs
+WangLi3SMil_F
+```
+
+## Understanding Commutative Noise
+
+### Commutativity Condition
+
+Noise terms g₁, g₂, ..., gₘ are **commutative** if:
+
+```
+[gᵢ, gⱼ] = gᵢ(∂gⱼ/∂x) - gⱼ(∂gᵢ/∂x) = 0
+```
+
+for all pairs (i,j).
+
+### Examples of Commutative Noise:
+
+ 1. **Additive noise**: g(u,t) = σ(t) (independent of u)
+ 2. **Scalar multiplicative**: g(u,t) = σ(t)u (same u dependence)
+ 3. **Diagonal with same function**: gᵢ(u,t) = σᵢ(t)h(u)
+
+### Examples of Non-commutative Noise:
+
+ 1. **Different multiplicative terms**: g₁ = σ₁u₁, g₂ = σ₂u₂
+ 2. **Cross-coupling**: g₁ = σ₁₁u₁ + σ₁₂u₂, g₂ = σ₂₁u₁ + σ₂₂u₂
+
+## Method Selection Guide
+
+### For Commutative Noise:
+
+ 1. **RKMilCommute** - Adaptive, excellent general choice
+ 2. **WangLi3SMil methods** - Fixed step, when dt is predetermined
+
+### For Non-commutative Noise:
+
+ 1. **RKMilGeneral** - Handles general case with Lévy area
+ 2. **Fall back to SRI/SRA methods** - More robust for complex noise
+
+### For Uncertain Commutativity:
+
+  - Test with RKMilCommute first
+  - If results seem incorrect, switch to RKMilGeneral or SRI methods
+
+## Computational Advantages
+
+**Commutative case:**
+
+  - No Lévy area computation needed
+  - Simpler stochastic integrals
+  - Higher efficiency per step
+  - Better scalability to high dimensions
+
+**Non-commutative case:**
+
+  - Requires Lévy area approximation
+  - More expensive per step
+  - Uses specialized algorithms (LevyArea.jl integration)
+
+## Performance Tips
+
+ 1. **Verify commutativity** before using specialized methods
+ 2. **Use appropriate interpretation** (Itô vs Stratonovich)
+ 3. **Consider problem dimension** - benefits increase with system size
+ 4. **Test accuracy** - commutative methods may be more sensitive
+
+## Iterated Integrals and Lévy Area
+
+### For Commutative Noise (IICommutative):
+
+Only simple stochastic integrals ∫₀ᵗ dWₛ are needed.
+
+### For Non-commutative Noise (IILevyArea):
+
+Requires double integrals ∫₀ᵗ ∫₀ˢ dWᵤdWₛ (Lévy area).
+
+RKMilGeneral automatically chooses optimal Lévy area algorithms via LevyArea.jl.
+
+## References
+
+  - Kloeden, P.E., Platen, E., "Numerical Solution of Stochastic Differential Equations"
+  - Wang and Li, "Three-stage stochastic Runge-Kutta methods for stochastic differential equations"
+  - Kastner, F. and Rößler, A., "LevyArea.jl" for Lévy area computation
diff --git a/lib/StochasticDiffEq/docs/src/nonstiff/high_weak_order.md b/lib/StochasticDiffEq/docs/src/nonstiff/high_weak_order.md
new file mode 100644
index 00000000000..1ff3871633e
--- /dev/null
+++ b/lib/StochasticDiffEq/docs/src/nonstiff/high_weak_order.md
@@ -0,0 +1,149 @@
+# High Weak Order Methods
+
+These methods are specifically designed for problems where weak convergence is more important than strong convergence. They are optimal for Monte Carlo simulations, computing expectations, moments, and other statistical properties of solutions.
+
+## Recommended High Weak Order Methods
+
+### DRI1 - Debrabant-Rößler Method (Weak Order 2)
+
+```@docs
+DRI1
+```
+
+### DRI1NM - Debrabant-Rößler for Non-mixing Diagonal Problems
+
+```@docs
+DRI1NM
+```
+
+## Other Weak Order 2 Methods
+
+### RI1, RI3, RI5, RI6 - Rößler Methods
+
+```@docs
+RI1
+```
+
+```@docs
+RI3
+```
+
+```@docs
+RI5
+```
+
+```@docs
+RI6
+```
+
+### RDI Methods - Alternative Weak Order 2
+
+```@docs
+RDI2WM
+```
+
+```@docs
+RDI3WM
+```
+
+```@docs
+RDI4WM
+```
+
+### W2Ito1 - Efficient Weak Order 2
+
+```@docs
+W2Ito1
+```
+
+## Fixed Step Methods
+
+### PL1WM, PL1WMA - Platen Methods
+
+```@docs
+PL1WM
+```
+
+```@docs
+PL1WMA
+```
+
+### Stratonovich Methods
+
+#### RS1, RS2 - Rößler Stratonovich Methods
+
+```@docs
+RS1
+```
+
+```@docs
+RS2
+```
+
+#### NON, NON2 - Non-commutative Stratonovich
+
+```@docs
+NON
+```
+
+```@docs
+NON2
+```
+
+#### COM - Commutative Stratonovich
+
+```@docs
+COM
+```
+
+## Specialized Methods
+
+### SIEA, SMEA, SIEB, SMEB - Tocino-Vigo-Aguiar Methods
+
+```@docs
+SIEA
+```
+
+```@docs
+SMEA
+```
+
+```@docs
+SIEB
+```
+
+```@docs
+SMEB
+```
+
+## Weak vs Strong Convergence
+
+**Strong Convergence**: Measures pathwise error E[|X(T) - X_h(T)|^p]
+**Weak Convergence**: Measures error in expectations E[f(X(T))] - E[f(X_h(T))]
+
+### When to Use Weak Order Methods:
+
+  - Monte Carlo simulations
+  - Computing expectations and moments
+  - Statistical analysis of SDEs
+  - When pathwise accuracy is not critical
+  - Large ensemble simulations
+
+### Advantages:
+
+  - Often more efficient for statistical quantities
+  - Can use larger time steps while maintaining weak accuracy
+  - Optimized error constants for better practical performance
+
+## Method Selection Guide
+
+ 1. **General purpose weak order 2**: DRI1
+ 2. **Non-mixing diagonal**: DRI1NM
+ 3. **Fixed step**: PL1WM, RS1/RS2
+ 4. **Stratonovich**: RS1, RS2, NON, NON2
+ 5. **Specialized applications**: RI methods, RDI methods
+
+## References
+
+  - Debrabant, K. and Rößler A., "Families of efficient second order Runge–Kutta methods for the weak approximation of Itô stochastic differential equations"
+  - Rößler A., "Second Order Runge–Kutta Methods for Itô Stochastic Differential Equations"
diff --git a/lib/StochasticDiffEq/docs/src/nonstiff/sra_sri_methods.md b/lib/StochasticDiffEq/docs/src/nonstiff/sra_sri_methods.md
new file mode 100644
index 00000000000..194e768c9f0
--- /dev/null
+++ b/lib/StochasticDiffEq/docs/src/nonstiff/sra_sri_methods.md
@@ -0,0 +1,128 @@
+# SRA/SRI Methods - Stochastic Runge-Kutta
+
+The SRA (Stochastic Runge-Kutta for Additive noise) and SRI (Stochastic Runge-Kutta for Itô) methods provide high-order adaptive solvers for different noise structures. These are among the most effective methods for their respective problem classes.
+
+## Recommended Methods
+
+### SOSRI - Stability-Optimized SRI (Recommended)
+
+```@docs
+SOSRI
+```
+
+### SOSRA - Stability-Optimized SRA (Optimal for Additive Noise)
+
+```@docs
+SOSRA
+```
+
+## Alternative SRI Methods
+
+### SRIW1 - SRI Weak Order 2
+
+```@docs
+SRIW1
+```
+
+### SRIW2 - SRI Weak Order 3
+
+```@docs
+SRIW2
+```
+
+### SOSRI2 - Alternative Stability-Optimized SRI
+
+```@docs
+SOSRI2
+```
+
+## Alternative SRA Methods
+
+### SRA1 - Original SRA Method
+
+```@docs
+SRA1
+```
+
+### SRA2 - SRA Method Version 2
+
+```@docs
+SRA2
+```
+
+### SRA3 - SRA Method with Weak Order 3
+
+```@docs
+SRA3
+```
+
+### SOSRA2 - Alternative Stability-Optimized SRA
+
+```@docs
+SOSRA2
+```
+
+## Configurable Methods
+
+### SRA - Configurable SRA with Custom Tableaux
+
+```@docs
+SRA
+```
+
+### SRI - Configurable SRI with Custom Tableaux
+
+```@docs
+SRI
+```
+
+## Method Selection Guide
+
+### For Diagonal/Scalar Noise:
+
+ 1. **First choice: SOSRI** - Best overall performance and stability
+ 2. **Alternative: SRIW1** - Standard SRI method
+ 3. **High weak order: SRIW2** - When weak order 3 is needed
+
+### For Additive Noise:
+
+ 1. **First choice: SOSRA** - Optimal for additive noise structure
+ 2. **Alternative: SRA1** - Standard SRA method
+ 3. **High weak order: SRA3** - When weak order 3 is needed
+
+### Performance Characteristics:
+
+  - **SOSRI/SOSRA**: Stability-optimized, robust to high tolerances
+  - **SRIWx/SRAx**: Standard methods with proven theoretical properties
+  - **SRA/SRI**: Allow custom tableaux for specialized applications
+
+## Theoretical Foundation
+
+The SRA and SRI methods are based on stochastic Runge-Kutta theory:
+
+**SRA Methods** exploit the additive noise structure:
+
+```
+du = f(u,t)dt + σ(t)dW
+```
+
+Where the diffusion σ doesn't depend on the solution u.
+
+**SRI Methods** handle the general diagonal case:
+
+```
+du = f(u,t)dt + g(u,t)dW
+```
+
+Where each component has independent noise.
+
+Both method families achieve:
+
+  - Strong order 1.5 convergence
+  - Weak order 2.0 or higher
+  - Adaptive time stepping with embedded error estimation
+  - A-stable or L-stable properties (for optimized versions)
+
+## References
+
+  - Rößler A., "Runge–Kutta Methods for the Strong Approximation of Solutions of Stochastic Differential Equations", SIAM J. Numer. Anal., 48 (3), pp. 922–952
diff --git a/lib/StochasticDiffEq/docs/src/stiff/implicit_methods.md b/lib/StochasticDiffEq/docs/src/stiff/implicit_methods.md
new file mode 100644
index 00000000000..1c3c31ab7ab
--- /dev/null
+++ b/lib/StochasticDiffEq/docs/src/stiff/implicit_methods.md
@@ -0,0 +1,108 @@
+# Implicit Methods for Stiff SDEs
+
+When SDEs have stiff drift terms, explicit methods may require impractically small time steps for stability. Implicit methods treat the drift term implicitly while keeping the diffusion explicit, providing excellent stability properties.
+
+## Recommended Stiff Methods
+
+### SKenCarp - Stochastic KenCarp (Highly Recommended for Stiff Additive Noise)
+
+```@docs
+SKenCarp
+```
+
+## Basic Implicit Methods
+
+### ImplicitEM - Implicit Euler-Maruyama
+
+```@docs
+ImplicitEM
+```
+
+### ImplicitEulerHeun - Implicit Euler-Heun (Stratonovich)
+
+```@docs
+ImplicitEulerHeun
+```
+
+### ImplicitRKMil - Implicit Runge-Kutta Milstein
+
+```@docs
+ImplicitRKMil
+```
+
+## Split-Step Implicit Methods
+
+### ISSEM - Implicit Split-Step Euler-Maruyama
+
+```@docs
+ISSEM
+```
+
+### ISSEulerHeun - Implicit Split-Step Euler-Heun
+
+```@docs
+ISSEulerHeun
+```
+
+## Method Selection Guide
+
+### Problem Classification:
+
+ 1. **Mildly stiff drift**: ImplicitEM, ImplicitRKMil
+ 2. **Stiff additive noise**: SKenCarp (highly recommended)
+ 3. **Fully stiff (including diffusion)**: ISSEM, ISSEulerHeun
+ 4. **Stratonovich problems**: ImplicitEulerHeun, ISSEulerHeun
+
+### Performance Ranking for Stiff Problems:
+
+ 1. **SKenCarp** - Best for stiff problems with additive noise
+ 2. **ISSEM/ISSEulerHeun** - For fully stiff problems
+ 3. **ImplicitRKMil** - Higher order for mildly stiff problems
+ 4. **ImplicitEM** - Robust fallback option
+
+## Understanding Stiffness in SDEs
+
+**Drift Stiffness**: Large negative eigenvalues in the drift term f(u,t) requiring small time steps for explicit stability.
+
+**Diffusion Stiffness**: Large coefficients in the diffusion term g(u,t) causing stability issues.
+
+**Detection**: If explicit methods require very small dt or produce unstable solutions, try implicit methods.
+
+## Configuration Options
+
+All implicit methods share common configuration parameters:
+
+```julia
+# Linear solver options
+ImplicitEM(linsolve = KrylovJL_GMRES())
+
+# Nonlinear solver options  
+ImplicitEM(nlsolve = NLNewton(max_iter = 20))
+
+# Jacobian computation
+ImplicitEM(autodiff = true, chunk_size = 4)
+
+# Theta method parameter
+ImplicitEM(theta = 0.5)  # Trapezoidal rule
+```
+
+## Symplectic Methods
+
+For Hamiltonian SDEs, use symplectic variants:
+
+```julia
+SImplicitMidpoint()  # Symplectic implicit midpoint
+STrapezoid()        # Symplectic trapezoidal rule
+```
+
+## Performance Tips
+
+ 1. **Jacobian**: Provide analytical Jacobian when possible
+ 2. **Linear solver**: Choose appropriate solver for problem structure
+ 3. **Preconditioning**: Use preconditioners for large systems
+ 4. **Theta parameter**: θ=0.5 often provides good accuracy/stability balance
+
+## References
+
+  - Standard implicit ODE methods adapted to SDEs
+  - Milstein, G.N., "Numerical Integration of Stochastic Differential Equations"
diff --git a/lib/StochasticDiffEq/docs/src/stiff/split_step_methods.md b/lib/StochasticDiffEq/docs/src/stiff/split_step_methods.md
new file mode 100644
index 00000000000..df3ee1de8a9
--- /dev/null
+++ b/lib/StochasticDiffEq/docs/src/stiff/split_step_methods.md
@@ -0,0 +1,87 @@
+# Split-Step Methods for Fully Stiff Problems
+
+When both drift and diffusion terms are stiff, split-step methods treat both parts implicitly. These methods can handle the most challenging stiff SDE problems.
+
+## Split-Step Implicit Methods
+
+### ISSEM - Implicit Split-Step Euler-Maruyama
+
+```@docs
+ISSEM
+```
+
+### ISSEulerHeun - Implicit Split-Step Euler-Heun
+
+```@docs
+ISSEulerHeun
+```
+
+## When to Use Split-Step Methods
+
+**Use ISSEM when:**
+
+  - Both drift and diffusion are stiff
+  - Explicit methods fail even with small time steps
+  - Standard implicit methods (ImplicitEM) are insufficient
+  - Working with Itô interpretation
+
+**Use ISSEulerHeun when:**
+
+  - Both drift and diffusion are stiff
+  - Working with Stratonovich interpretation
+  - Need implicit treatment of diffusion term
+
+## Understanding Full Stiffness
+
+**Drift Stiffness**: Large negative eigenvalues in f(u,t)
+**Diffusion Stiffness**: Large coefficients in g(u,t) causing instability
+
+**Detection Signs:**
+
+  - ImplicitEM still requires very small time steps
+  - Solutions become unstable despite implicit drift treatment
+  - Large diffusion coefficients cause numerical artifacts
+
+## Algorithm Structure
+
+Split-step methods solve:
+
+```
+du = f(u,t)dt + g(u,t)dW
+```
+
+By treating both terms implicitly through operator splitting or fully implicit schemes.
+
+## Performance Considerations
+
+  - More expensive per step than drift-only implicit methods
+  - May require smaller time steps than expected
+  - Jacobian computations for both drift and diffusion
+  - Nonlinear solve complexity increases
+
+## Configuration
+
+Same options as other implicit methods:
+
+```julia
+ISSEM(
+    linsolve = KrylovJL_GMRES(),
+    nlsolve = NLNewton(),
+    theta = 1.0,
+    autodiff = true
+)
+```
+
+## Alternative Approaches
+
+If split-step methods are too expensive:
+
+ 1. Try stabilized explicit methods (SROCK family)
+ 2. Consider method of lines for PDE problems
+ 3. Use shorter time intervals with restarts
+ 4. Reformulate problem to reduce stiffness
+
+## References
+
+  - Implicit-explicit splitting schemes for SDEs
+  - Stochastic operator splitting methods
diff --git a/lib/StochasticDiffEq/docs/src/stiff/stabilized_methods.md b/lib/StochasticDiffEq/docs/src/stiff/stabilized_methods.md
new file mode 100644
index 00000000000..5b4af39cb7d
--- /dev/null
+++ b/lib/StochasticDiffEq/docs/src/stiff/stabilized_methods.md
@@ -0,0 +1,138 @@
+# Stabilized Methods (SROCK Family)
+
+Stabilized Runge-Kutta Chebyshev (SROCK) methods provide stability for mildly stiff problems through extended stability regions rather than implicit treatment. These methods are particularly effective for parabolic PDEs discretized by method of lines.
+
+## SROCK Methods
+
+### SROCK1 - First Order Stabilized Method
+
+```@docs
+SROCK1
+```
+
+### SROCK2, KomBurSROCK2, SROCKC2 - Second Order Methods
+
+```@docs
+SROCK2
+```
+
+```@docs
+KomBurSROCK2
+```
+
+```@docs
+SROCKC2
+```
+
+### SROCKEM - Stabilized Euler-Maruyama
+
+```@docs
+SROCKEM
+```
+
+### SKSROCK - Stabilized Method with Post-Processing
+
+```@docs
+SKSROCK
+```
+
+### TangXiaoSROCK2 - Alternative Second Order Method
+
+```@docs
+TangXiaoSROCK2
+```
+
+## When to Use Stabilized Methods
+
+**Ideal for:**
+
+  - Parabolic PDEs discretized by method of lines
+  - Problems with moderate stiffness (not extremely stiff)
+  - Large systems where implicit methods are expensive
+  - When eigenvalue spectrum is primarily negative real
+
+**Advantages over implicit methods:**
+
+  - No linear system solves required
+  - Better for large systems
+  - Parallelizable
+  - No Jacobian computation needed
+
+**Disadvantages:**
+
+  - Limited to moderate stiffness
+  - May require eigenvalue estimation
+  - Not effective for highly oscillatory problems
+
+## Stability Regions
+
+SROCK methods extend stability along the negative real axis:
+
+  - **Standard explicit**: Stability region ~ [-2, 0]
+  - **SROCK methods**: Stability region ~ [-s², 0] where s is the number of stages
+
+The number of stages s is chosen based on estimated eigenvalues.
+
+## Eigenvalue Estimation
+
+Most SROCK methods accept an `eigen_est` parameter:
+
+```julia
+# Automatic estimation (default)
+SROCK1()
+
+# Manual estimation
+SROCK1(eigen_est = -100.0)  # Largest eigenvalue magnitude
+
+# Custom estimation function
+SROCK1(eigen_est = my_estimator)
+```
+
+## Method Selection Guide
+
+ 1. **SROCK1**: Basic first-order method, most robust
+ 2. **SROCK2**: Higher accuracy, good general choice
+ 3. **SROCKEM**: When Euler-Maruyama structure is preferred
+ 4. **SKSROCK**: Advanced features, post-processing options
+ 5. **SROCKC2**: Conservative second-order variant
+
+## Problem Suitability
+
+**Well-suited:**
+
+  - Reaction-diffusion equations
+  - Heat equations with stochastic terms
+  - Large sparse systems
+  - Method of lines discretizations
+
+**Not well-suited:**
+
+  - Highly stiff problems (use implicit methods)
+  - Problems with complex eigenvalue spectra
+  - Small dense systems (overhead not justified)
+
+## Configuration Tips
+
+```julia
+# For PDE problems
+SROCK2(eigen_est = estimate_spectral_radius(A))
+
+# For uncertain problems, start conservatively
+SROCK1()  # Most robust
+
+# For higher accuracy
+SROCK2()  # Good balance
+```
+
+## Performance Considerations
+
+  - Stage count increases with stiffness
+  - Eigenvalue estimation cost
+  - Memory requirements for internal stages
+  - Better scalability than implicit methods
+
+## References
+
+  - Chebyshev methods for parabolic problems
+  - ROCK methods for stiff ODEs
+  - Stabilized explicit methods for PDEs
diff --git a/lib/StochasticDiffEq/docs/src/usage.md b/lib/StochasticDiffEq/docs/src/usage.md
new file mode 100644
index 00000000000..ae40ecd8935
--- /dev/null
+++ b/lib/StochasticDiffEq/docs/src/usage.md
@@ -0,0 +1,239 @@
+# Usage Guide
+
+This page provides guidance on using StochasticDiffEq.jl effectively.
+
+## Basic Usage
+
+### Problem Definition
+
+StochasticDiffEq.jl uses the standard DifferentialEquations.jl problem interface:
+
+```julia
+using StochasticDiffEq
+
+# For scalar problems
+function f(u, p, t)  # drift
+    return μ * u
+end
+
+function g(u, p, t)  # diffusion  
+    return σ * u
+end
+
+# For in-place systems
+function f!(du, u, p, t)
+    du[1] = μ * u[1]
+    du[2] = -ν * u[2]
+end
+
+function g!(du, u, p, t)
+    du[1] = σ₁ * u[1]
+    du[2] = σ₂ * u[2]
+end
+
+# Create problem
+prob = SDEProblem(f, g, u0, tspan)
+```
+
+### Solver Selection
+
+Choose solvers based on your problem characteristics:
+
+```julia
+# Default - good for most problems
+sol = solve(prob)
+
+# Specify solver explicitly
+sol = solve(prob, SOSRI())          # Recommended for diagonal noise
+sol = solve(prob, SOSRA())          # Optimal for additive noise  
+sol = solve(prob, SKenCarp())       # For stiff problems
+sol = solve(prob, EM())             # For maximum efficiency
+```
+
+### Algorithm Parameters
+
+Most solvers accept parameters for customization:
+
+```julia
+# Euler-Maruyama with step splitting
+sol = solve(prob, EM(split = true))
+
+# RKMilCommute with Stratonovich interpretation
+sol = solve(prob, RKMilCommute(interpretation = :Stratonovich))
+
+# Implicit methods with solver options
+sol = solve(prob, SKenCarp(linsolve = KrylovJL_GMRES()))
+```
+
+## Tolerances and Adaptive Stepping
+
+Set absolute and relative tolerances:
+
+```julia
+sol = solve(prob, SOSRI(), abstol = 1e-6, reltol = 1e-3)
+```
+
+For fixed time stepping:
+
+```julia
+sol = solve(prob, EM(), dt = 0.01, adaptive = false)
+```
+
+## Noise Types
+
+### Diagonal Noise
+
+Most common case - each component has independent noise:
+
+```julia
+function g!(du, u, p, t)
+    du[1] = σ₁ * u[1]
+    du[2] = σ₂ * u[2]
+end
+```
+
+### Scalar Noise
+
+Single noise source affects all components:
+
+```julia
+function g!(du, u, p, t)
+    du[1] = σ * u[1]
+    du[2] = σ * u[2]
+end
+```
+
+### Non-diagonal Noise
+
+Multiple noise sources with cross-terms:
+
+```julia
+function g!(du, u, p, t)
+    du[1] = σ₁₁ * u[1] + σ₁₂ * u[2]
+    du[2] = σ₂₁ * u[1] + σ₂₂ * u[2]
+end
+```
+
+### Additive Noise
+
+Noise independent of solution:
+
+```julia
+function g!(du, u, p, t)
+    du[1] = σ₁
+    du[2] = σ₂
+end
+```
+
+## Itô vs Stratonovich
+
+Specify interpretation when creating problems or choosing solvers:
+
+```julia
+# Itô interpretation (default)
+prob = SDEProblem(f!, g!, u0, tspan, interpretation = :Ito)
+
+# Stratonovich interpretation  
+prob = SDEProblem(f!, g!, u0, tspan, interpretation = :Stratonovich)
+
+# Or at solver level
+sol = solve(prob, RKMil(interpretation = :Stratonovich))
+```
+
+## RNG Control
+
+Each `solve` or `init` call needs a random number generator (RNG) for
+constructing noise processes. The RNG is resolved from the keyword arguments
+in the following priority order:
+
+1. **`rng` provided** — use the given `AbstractRNG` directly. Any `seed` kwarg
+   is ignored. If a `TaskLocalRNG` is passed (i.e. `Random.default_rng()`), it
+   is converted to a concrete `Xoshiro` seeded from one draw of the task-local
+   stream, so the integrator never shares the global random stream.
+2. **`seed` provided (nonzero)** — construct `Xoshiro(seed)`.
+3. **Problem seed** (`prob.seed != 0`) — construct `Xoshiro(prob.seed)`.
+4. **Neither** — generate a random seed and construct `Xoshiro` from it.
+
+### Examples
+
+```julia
+using Random
+
+# Reproducible results with an explicit RNG
+rng = Xoshiro(42)
+sol = solve(prob, EM(); dt = 0.01, rng)
+
+# Same seed produces identical trajectories
+rng2 = Xoshiro(42)
+sol2 = solve(prob, EM(); dt = 0.01, rng = rng2)
+sol.u == sol2.u  # true
+
+# The older seed keyword still works (constructs Xoshiro(seed) internally,
+# so reproducibility depends on the internal RNG type; prefer `rng` for
+# guaranteed reproducibility across library versions)
+sol = solve(prob, EM(); dt = 0.01, seed = UInt64(42))
+```
+
+### RNG ownership
+
+The `rng` keyword controls **framework-constructed** randomness (noise processes
+created internally by the solver, and the integrator's own RNG). If you supply
+your own noise process via `SDEProblem(...; noise = my_W)`, that noise object's
+internal RNG remains under your control and is **not** modified by the `rng`
+keyword or by `set_rng!`.
+
+### Integrator RNG interface
+
+The integrator implements the SciMLBase RNG interface:
+
+```julia
+integ = init(prob, EM(); dt = 0.01, rng = Xoshiro(42))
+
+SciMLBase.has_rng(integ)       # true
+SciMLBase.get_rng(integ)       # the Xoshiro RNG
+SciMLBase.set_rng!(integ, Xoshiro(99))  # replace with same-type RNG
+```
+
+`set_rng!` requires the new RNG to be the same concrete type as the current one.
+For framework-constructed noise, it also syncs the noise process RNG. For
+user-provided noise, only `integrator.rng` is updated.
+
+The `reinit!` function also accepts an `rng` keyword:
+
+```julia
+reinit!(integ, u0; rng = Xoshiro(99))
+```
+
+## Performance Tips
+
+ 1. **Use appropriate solvers**: Match solver to problem type
+ 2. **In-place functions**: Use `f!(du,u,p,t)` for better performance
+ 3. **Tolerances**: Don't make tolerances unnecessarily strict
+ 4. **Static arrays**: Use `StaticArrays.jl` for small systems
+ 5. **GPU**: Use `CuArrays.jl` for large problems
+
+## Common Pitfalls
+
+ 1. **Wrong noise type**: Ensure solver supports your noise structure
+ 2. **Stiffness**: Use appropriate stiff solvers for stiff problems
+ 3. **Commuting noise**: Use specialized solvers for better efficiency
+ 4. **High dimensions**: Consider weak convergence methods for Monte Carlo
+
+## Integration with DifferentialEquations.jl
+
+StochasticDiffEq.jl integrates with the broader ecosystem:
+
+```julia
+using DifferentialEquations
+
+# Callbacks
+condition(u, t, integrator) = u[1] - 0.5
+affect!(integrator) = terminate!(integrator)
+cb = ContinuousCallback(condition, affect!)
+
+sol = solve(prob, SOSRI(), callback = cb)
+
+# Ensemble simulations
+monte_prob = EnsembleProblem(prob)
+sim = solve(monte_prob, SOSRI(), EnsembleThreads(), trajectories = 1000)
+```