{"id":106719,"date":"2026-09-18T06:40:32","date_gmt":"2026-09-17T22:40:32","guid":{"rendered":"https:\/\/www.wsisp.com\/helps\/106719.html"},"modified":"2026-09-18T06:40:32","modified_gmt":"2026-09-17T22:40:32","slug":"c%e4%bb%a3%e7%a0%81%e5%ae%9e%e7%8e%b0matlab%e4%b8%ad%e7%9a%84ode45%e5%87%bd%e6%95%b0%e5%8a%9f%e8%83%bd","status":"publish","type":"post","link":"https:\/\/www.wsisp.com\/helps\/106719.html","title":{"rendered":"C++\u4ee3\u7801\u5b9e\u73b0MATLAB\u4e2d\u7684ode45\u51fd\u6570\u529f\u80fd"},"content":{"rendered":"

#<\/span>include<\/span> <stddef.h><\/span><\/span><\/p>\n

#<\/span>ifdef<\/span> __cplusplus<\/span><\/span>
\nextern<\/span> "C"<\/span> {<\/span>
\n#<\/span>endif<\/span><\/span><\/p>\n

\/\/ C callback for derivatives: fills dydt of length n given y of length n at time t.<\/span>
\ntypedef<\/span> void<\/span> (<\/span>*<\/span>ode45_deriv_cb)<\/span>(<\/span>double<\/span> t,<\/span> const<\/span> double<\/span>*<\/span> y,<\/span> double<\/span>*<\/span> dydt,<\/span> size_t n,<\/span> void<\/span>*<\/span> ctx)<\/span>;<\/span><\/p>\n

\/\/ Integrate from t0 to t1 using Dormand-Prince 5(4) with adaptive steps<\/span><\/p>\n

\/\/ Integrate from t0 to t1 using Dormand-Prince 5(4) with adaptive steps.<\/span>
\n\/\/ Inputs:<\/span>
\n\/\/ – t0, y0[n]: initial time\/state<\/span>
\n\/\/ – t1: target time<\/span>
\n\/\/ – h_init: initial step size (non-zero)<\/span>
\n\/\/ – atol, rtol: tolerances<\/span>
\n\/\/ – normControl: 0 (component-wise max) or 1 (global 2-norm)<\/span>
\n\/\/ – max_steps: safety cap<\/span>
\n\/\/ – cb, ctx: derivative callback + user context<\/span>
\n\/\/ Output:<\/span>
\n\/\/ – y_out[n]: state at t1<\/span>
\n\/\/ Returns 0 on success, non-zero on failure.<\/span>
\nint<\/span> ode45_integrate_final<\/span>(<\/span>double<\/span> t0,<\/span>
\n const<\/span> double<\/span>*<\/span> y0,<\/span>
\n size_t n,<\/span>
\n double<\/span> t1,<\/span>
\n double<\/span> h_init,<\/span>
\n double<\/span> atol,<\/span>
\n double<\/span> rtol,<\/span>
\n int<\/span> normControl,<\/span>
\n int<\/span> max_steps,<\/span>
\n ode45_deriv_cb cb,<\/span>
\n void<\/span>*<\/span> ctx,<\/span>
\n double<\/span>*<\/span> y_out)<\/span>;<\/span><\/p>\n

\/\/ Integrate and return samples at explicit times.<\/span>
\n\/\/ Inputs:<\/span>
\n\/\/ – t0, y0[n]<\/span>
\n\/\/ – times[m]: strictly monotone in desired direction (first equals t0 recommended)<\/span>
\n\/\/ – Remaining parameters like above.<\/span>
\n\/\/ Output:<\/span>
\n\/\/ – y_out[m*n]: concatenated row-major blocks of length n for each time.<\/span>
\n\/\/ Returns 0 on success, non-zero on failure.<\/span>
\nint<\/span> ode45_integrate_times<\/span>(<\/span>double<\/span> t0,<\/span>
\n const<\/span> double<\/span>*<\/span> y0,<\/span>
\n size_t n,<\/span>
\n const<\/span> double<\/span>*<\/span> times,<\/span>
\n size_t m,<\/span>
\n double<\/span> h_init,<\/span>
\n double<\/span> atol,<\/span>
\n double<\/span> rtol,<\/span>
\n int<\/span> normControl,<\/span>
\n int<\/span> max_steps,<\/span>
\n ode45_deriv_cb cb,<\/span>
\n void<\/span>*<\/span> ctx,<\/span>
\n double<\/span>*<\/span> y_out)<\/span>;<\/span><\/p>\n

#<\/span>ifdef<\/span> __cplusplus<\/span><\/span>
\n}<\/span>
\n#<\/span>endif<\/span><\/span><\/p>\n

#<\/span>include<\/span> <cmath><\/span><\/span>
\n#<\/span>include<\/span> <cstdlib><\/span><\/span>
\n#<\/span>include<\/span> <functional><\/span><\/span>
\n#<\/span>include<\/span> <iomanip><\/span><\/span>
\n#<\/span>include<\/span> <iostream><\/span><\/span>
\n#<\/span>include<\/span> <limits><\/span><\/span>
\n#<\/span>include<\/span> <stdexcept><\/span><\/span>
\n#<\/span>include<\/span> <array><\/span><\/span>
\n#<\/span>include<\/span> <utility><\/span><\/span>
\n#<\/span>include<\/span> <string><\/span><\/span>
\n#<\/span>include<\/span> <vector><\/span><\/span><\/p>\n

\/\/ Standalone Dormand-Prince 5(4) (DOPRI5) adaptive Runge-Kutta integrator.<\/span>
\n\/\/ No Boost dependencies.<\/span>
\nusing<\/span> State =<\/span> std::<\/span>vector<<\/span>double<\/span>><\/span>;<\/span>
\nusing<\/span> DerivFunc =<\/span> std::<\/span>function<<\/span>void<\/span>(<\/span>double<\/span> \/*t*\/<\/span>,<\/span> const<\/span> State&<\/span> \/*y*\/<\/span>,<\/span> State&<\/span> \/*dydt*\/<\/span>)<\/span>><\/span>;<\/span><\/p>\n

\/\/ Forward declaration for NormControl two-norm error<\/span>
\nstatic<\/span> double<\/span> two_norm_scaled<\/span>(<\/span>const<\/span> State&<\/span> err,<\/span> const<\/span> State&<\/span> y,<\/span> const<\/span> State&<\/span> y_next,<\/span> double<\/span> atol,<\/span> double<\/span> rtol)<\/span>;<\/span>
\nstatic<\/span> State zeros_like<\/span>(<\/span>const<\/span> State&<\/span> y)<\/span> {<\/span>
\n return<\/span> State<\/span>(<\/span>y.<\/span>size<\/span>(<\/span>)<\/span>,<\/span> 0.0<\/span>)<\/span>;<\/span>
\n}<\/span><\/p>\n

\/\/ Dense output for Dormand-Prince 5(4): Hermite quartic interpolation (Octave-like)<\/span>
\nstatic<\/span> State dopri5_dense_output<\/span>(<\/span>double<\/span> dt,<\/span> const<\/span> State&<\/span> x0,<\/span> const<\/span> State&<\/span> x1,<\/span> const<\/span> std::<\/span>array<<\/span>State,<\/span>7<\/span>><\/span>&<\/span> K,<\/span> double<\/span> s)<\/span> {<\/span>
\n \/\/ s in [0,1], dt = h, K[0]=k1 at t, K[6]=k7 at t+dt<\/span>
\n const<\/span> size_t n =<\/span> x0.<\/span>size<\/span>(<\/span>)<\/span>;<\/span>
\n State out<\/span>(<\/span>n)<\/span>;<\/span>
\n if<\/span> (<\/span>n ==<\/span> 0<\/span>)<\/span> return<\/span> out;<\/span><\/p>\n

\/\/ Coefficients for u_half = y(t+dt\/2)<\/span>
\n constexpr<\/span> double<\/span> c_half[<\/span>7<\/span>]<\/span> =<\/span> {<\/span>
\n 6025192743.0<\/span>\/<\/span>30085553152.0<\/span>,<\/span>
\n 0.0<\/span>,<\/span>
\n 51252292925.0<\/span>\/<\/span>65400821598.0<\/span>,<\/span>
\n –<\/span>2691868925.0<\/span>\/<\/span>45128329728.0<\/span>,<\/span>
\n 187940372067.0<\/span>\/<\/span>1594534317056.0<\/span>,<\/span>
\n –<\/span>1776094331.0<\/span>\/<\/span>19743644256.0<\/span>,<\/span>
\n 11237099.0<\/span>\/<\/span>235043384.0<\/span>,<\/span>
\n }<\/span>;<\/span><\/p>\n

State u_half<\/span>(<\/span>n,<\/span> 0.0<\/span>)<\/span>;<\/span>
\n for<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> n;<\/span> ++<\/span>i)<\/span> {<\/span>
\n double<\/span> acc =<\/span> 0.0<\/span>;<\/span>
\n for<\/span> (<\/span>int<\/span> k =<\/span> 0<\/span>;<\/span> k <<\/span> 7<\/span>;<\/span> ++<\/span>k)<\/span> acc +=<\/span> c_half[<\/span>k]<\/span> *<\/span> K[<\/span>k]<\/span>[<\/span>i]<\/span>;<\/span>
\n u_half[<\/span>i]<\/span> =<\/span> x0[<\/span>i]<\/span> +<\/span> 0.5<\/span> *<\/span> dt *<\/span> acc;<\/span>
\n }<\/span><\/p>\n

\/\/ Hermite quartic basis polynomials<\/span>
\n const<\/span> double<\/span> s2 =<\/span> s*<\/span>s;<\/span>
\n const<\/span> double<\/span> s3 =<\/span> s2*<\/span>s;<\/span>
\n const<\/span> double<\/span> s4 =<\/span> s3*<\/span>s;<\/span><\/p>\n

const<\/span> double<\/span> H0 =<\/span> 1.0<\/span> –<\/span> 11.0<\/span>*<\/span>s2 +<\/span> 18.0<\/span>*<\/span>s3 –<\/span> 8.0<\/span>*<\/span>s4;<\/span>
\n const<\/span> double<\/span> H1 =<\/span> s –<\/span> 4.0<\/span>*<\/span>s2 +<\/span> 5.0<\/span>*<\/span>s3 –<\/span> 2.0<\/span>*<\/span>s4;<\/span>
\n const<\/span> double<\/span> H2 =<\/span> 16.0<\/span>*<\/span>s2 –<\/span> 32.0<\/span>*<\/span>s3 +<\/span> 16.0<\/span>*<\/span>s4;<\/span>
\n const<\/span> double<\/span> H3 =<\/span> –<\/span>5.0<\/span>*<\/span>s2 +<\/span> 14.0<\/span>*<\/span>s3 –<\/span> 8.0<\/span>*<\/span>s4;<\/span>
\n const<\/span> double<\/span> H4 =<\/span> s2 –<\/span> 3.0<\/span>*<\/span>s3 +<\/span> 2.0<\/span>*<\/span>s4;<\/span><\/p>\n

for<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> n;<\/span> ++<\/span>i)<\/span> {<\/span>
\n out[<\/span>i]<\/span> =<\/span> H0 *<\/span> x0[<\/span>i]<\/span>
\n +<\/span> H1 *<\/span> (<\/span>dt *<\/span> K[<\/span>0<\/span>]<\/span>[<\/span>i]<\/span>)<\/span>
\n +<\/span> H2 *<\/span> u_half[<\/span>i]<\/span>
\n +<\/span> H3 *<\/span> x1[<\/span>i]<\/span>
\n +<\/span> H4 *<\/span> (<\/span>dt *<\/span> K[<\/span>6<\/span>]<\/span>[<\/span>i]<\/span>)<\/span>;<\/span>
\n }<\/span><\/p>\n

return<\/span> out;<\/span>
\n}<\/span><\/p>\n

static<\/span> double<\/span> inf_norm_scaled<\/span>(<\/span>const<\/span> State&<\/span> err,<\/span> const<\/span> State&<\/span> y,<\/span> const<\/span> State&<\/span> y_next,<\/span> double<\/span> atol,<\/span> double<\/span> rtol)<\/span> {<\/span>
\n \/\/ Matches Octave's AbsRel_norm with NormControl == false:<\/span>
\n \/\/ err_norm = max_i (|err_i| \/ scale_i), scale_i = max(AbsTol, RelTol * max(|y_i|, |y_next_i|))<\/span>
\n if<\/span> (<\/span>err.<\/span>size<\/span>(<\/span>)<\/span> !=<\/span> y.<\/span>size<\/span>(<\/span>)<\/span> ||<\/span> err.<\/span>size<\/span>(<\/span>)<\/span> !=<\/span> y_next.<\/span>size<\/span>(<\/span>)<\/span>)<\/span> {<\/span>
\n throw<\/span> std::<\/span>runtime_error<\/span>(<\/span>"State size mismatch"<\/span>)<\/span>;<\/span>
\n }<\/span>
\n if<\/span> (<\/span>err.<\/span>empty<\/span>(<\/span>)<\/span>)<\/span> {<\/span>
\n return<\/span> 0.0<\/span>;<\/span>
\n }<\/span>
\n double<\/span> max_ratio =<\/span> 0.0<\/span>;<\/span>
\n for<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> err.<\/span>size<\/span>(<\/span>)<\/span>;<\/span> ++<\/span>i)<\/span> {<\/span>
\n const<\/span> double<\/span> sc =<\/span> std::<\/span>max<\/span>(<\/span>atol,<\/span> rtol *<\/span> std::<\/span>max<\/span>(<\/span>std::<\/span>abs<\/span>(<\/span>y[<\/span>i]<\/span>)<\/span>,<\/span> std::<\/span>abs<\/span>(<\/span>y_next[<\/span>i]<\/span>)<\/span>)<\/span>)<\/span>;<\/span>
\n const<\/span> double<\/span> r =<\/span> (<\/span>sc ><\/span> 0.0<\/span>)<\/span> ?<\/span> (<\/span>std::<\/span>abs<\/span>(<\/span>err[<\/span>i]<\/span>)<\/span> \/<\/span> sc)<\/span> :<\/span> std::<\/span>abs<\/span>(<\/span>err[<\/span>i]<\/span>)<\/span>;<\/span>
\n if<\/span> (<\/span>r ><\/span> max_ratio)<\/span> max_ratio =<\/span> r;<\/span>
\n }<\/span>
\n return<\/span> max_ratio;<\/span>
\n}<\/span><\/p>\n

\/\/ Performs one DOPRI5 step. Computes y_next (5th order) and err (difference 5th-4th).<\/span>
\nstatic<\/span> void<\/span> runge_kutta_dopri5_step<\/span>(<\/span>const<\/span> DerivFunc&<\/span> f,<\/span> double<\/span> t,<\/span> const<\/span> State&<\/span> y,<\/span> double<\/span> h,<\/span> State&<\/span> y_next,<\/span>
\nState&<\/span> err,<\/span> std::<\/span>array<<\/span>State,<\/span>7<\/span>><\/span>&<\/span> K)<\/span> {<\/span>
\n const<\/span> size_t n =<\/span> y.<\/span>size<\/span>(<\/span>)<\/span>;<\/span>
\n if<\/span> (<\/span>n ==<\/span> 0<\/span>)<\/span> {<\/span>
\n y_next.<\/span>clear<\/span>(<\/span>)<\/span>;<\/span>
\n err.<\/span>clear<\/span>(<\/span>)<\/span>;<\/span>
\n return<\/span>;<\/span>
\n }<\/span><\/p>\n

State k1<\/span>(<\/span>n)<\/span>,<\/span> k2<\/span>(<\/span>n)<\/span>,<\/span> k3<\/span>(<\/span>n)<\/span>,<\/span> k4<\/span>(<\/span>n)<\/span>,<\/span> k5<\/span>(<\/span>n)<\/span>,<\/span> k6<\/span>(<\/span>n)<\/span>,<\/span> k7<\/span>(<\/span>n)<\/span>;<\/span>
\n State tmp<\/span>(<\/span>n)<\/span>;<\/span><\/p>\n

\/\/ Butcher tableau coefficients for Dormand-Prince 5(4)<\/span>
\n \/\/ c_i<\/span>
\n constexpr<\/span> double<\/span> c2 =<\/span> 1.0<\/span> \/<\/span> 5.0<\/span>;<\/span>
\n constexpr<\/span> double<\/span> c3 =<\/span> 3.0<\/span> \/<\/span> 10.0<\/span>;<\/span>
\n constexpr<\/span> double<\/span> c4 =<\/span> 4.0<\/span> \/<\/span> 5.0<\/span>;<\/span>
\n constexpr<\/span> double<\/span> c5 =<\/span> 8.0<\/span> \/<\/span> 9.0<\/span>;<\/span>
\n constexpr<\/span> double<\/span> c6 =<\/span> 1.0<\/span>;<\/span>
\n constexpr<\/span> double<\/span> c7 =<\/span> 1.0<\/span>;<\/span><\/p>\n

\/\/ a (lower triangular):<\/span>
\n constexpr<\/span> double<\/span> a21 =<\/span> 1.0<\/span> \/<\/span> 5.0<\/span>;<\/span><\/p>\n

constexpr<\/span> double<\/span> a31 =<\/span> 3.0<\/span> \/<\/span> 40.0<\/span>;<\/span>
\n constexpr<\/span> double<\/span> a32 =<\/span> 9.0<\/span> \/<\/span> 40.0<\/span>;<\/span><\/p>\n

constexpr<\/span> double<\/span> a41 =<\/span> 44.0<\/span> \/<\/span> 45.0<\/span>;<\/span>
\n constexpr<\/span> double<\/span> a42 =<\/span> –<\/span>56.0<\/span> \/<\/span> 15.0<\/span>;<\/span>
\n constexpr<\/span> double<\/span> a43 =<\/span> 32.0<\/span> \/<\/span> 9.0<\/span>;<\/span><\/p>\n

constexpr<\/span> double<\/span> a51 =<\/span> 19372.0<\/span> \/<\/span> 6561.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> a52 =<\/span> –<\/span>25360.0<\/span> \/<\/span> 2187.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> a53 =<\/span> 64448.0<\/span> \/<\/span> 6561.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> a54 =<\/span> –<\/span>212.0<\/span> \/<\/span> 729.0<\/span>;<\/span><\/p>\n

constexpr<\/span> double<\/span> a61 =<\/span> 9017.0<\/span> \/<\/span> 3168.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> a62 =<\/span> –<\/span>355.0<\/span> \/<\/span> 33.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> a63 =<\/span> 46732.0<\/span> \/<\/span> 5247.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> a64 =<\/span> 49.0<\/span> \/<\/span> 176.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> a65 =<\/span> –<\/span>5103.0<\/span> \/<\/span> 18656.0<\/span>;<\/span><\/p>\n

constexpr<\/span> double<\/span> a71 =<\/span> 35.0<\/span> \/<\/span> 384.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> a72 =<\/span> 0.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> a73 =<\/span> 500.0<\/span> \/<\/span> 1113.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> a74 =<\/span> 125.0<\/span> \/<\/span> 192.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> a75 =<\/span> –<\/span>2187.0<\/span> \/<\/span> 6784.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> a76 =<\/span> 11.0<\/span> \/<\/span> 84.0<\/span>;<\/span><\/p>\n

\/\/ b (5th order) equals row 7 coefficients (a71..a76) with b7=0.<\/span>
\nconstexpr<\/span> double<\/span> b1 =<\/span> a71;<\/span>
\nconstexpr<\/span> double<\/span> b2 =<\/span> a72;<\/span>
\nconstexpr<\/span> double<\/span> b3 =<\/span> a73;<\/span>
\nconstexpr<\/span> double<\/span> b4 =<\/span> a74;<\/span>
\nconstexpr<\/span> double<\/span> b5 =<\/span> a75;<\/span>
\nconstexpr<\/span> double<\/span> b6 =<\/span> a76;<\/span>
\nconstexpr<\/span> double<\/span> b7 =<\/span> 0.0<\/span>;<\/span><\/p>\n

\/\/ b_hat (4th order):<\/span>
\nconstexpr<\/span> double<\/span> bh1 =<\/span> 5179.0<\/span> \/<\/span> 57600.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> bh2 =<\/span> 0.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> bh3 =<\/span> 7571.0<\/span> \/<\/span> 16695.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> bh4 =<\/span> 393.0<\/span> \/<\/span> 640.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> bh5 =<\/span> –<\/span>92097.0<\/span> \/<\/span> 339200.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> bh6 =<\/span> 187.0<\/span> \/<\/span> 2100.0<\/span>;<\/span>
\nconstexpr<\/span> double<\/span> bh7 =<\/span> 1.0<\/span> \/<\/span> 40.0<\/span>;<\/span><\/p>\n

\/\/ k1<\/span>
\nf<\/span>(<\/span>t,<\/span> y,<\/span> k1)<\/span>;<\/span><\/p>\n

\/\/ k2<\/span>
\nfor<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> n;<\/span> ++<\/span>i)<\/span> tmp[<\/span>i]<\/span> =<\/span> y[<\/span>i]<\/span> +<\/span> h *<\/span> (<\/span>a21 *<\/span> k1[<\/span>i]<\/span>)<\/span>;<\/span>
\nf<\/span>(<\/span>t +<\/span> c2 *<\/span> h,<\/span> tmp,<\/span> k2)<\/span>;<\/span><\/p>\n

\/\/ k3<\/span>
\nfor<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> n;<\/span> ++<\/span>i)<\/span> tmp[<\/span>i]<\/span> =<\/span> y[<\/span>i]<\/span> +<\/span> h *<\/span> (<\/span>a31 *<\/span> k1[<\/span>i]<\/span> +<\/span> a32 *<\/span> k2[<\/span>i]<\/span>)<\/span>;<\/span>
\nf<\/span>(<\/span>t +<\/span> c3 *<\/span> h,<\/span> tmp,<\/span> k3)<\/span>;<\/span><\/p>\n

\/\/ k4<\/span>
\nfor<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> n;<\/span> ++<\/span>i)<\/span> tmp[<\/span>i]<\/span> =<\/span> y[<\/span>i]<\/span> +<\/span> h *<\/span> (<\/span>a41 *<\/span> k1[<\/span>i]<\/span> +<\/span> a42 *<\/span> k2[<\/span>i]<\/span> +<\/span> a43 *<\/span> k3[<\/span>i]<\/span>)<\/span>;<\/span>
\nf<\/span>(<\/span>t +<\/span> c4 *<\/span> h,<\/span> tmp,<\/span> k4)<\/span>;<\/span><\/p>\n

\/\/ k5<\/span>
\nfor<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> n;<\/span> ++<\/span>i)<\/span>
\n tmp[<\/span>i]<\/span> =<\/span> y[<\/span>i]<\/span> +<\/span> h *<\/span> (<\/span>a51 *<\/span> k1[<\/span>i]<\/span> +<\/span> a52 *<\/span> k2[<\/span>i]<\/span> +<\/span> a53 *<\/span> k3[<\/span>i]<\/span> +<\/span> a54 *<\/span> k4[<\/span>i]<\/span>)<\/span>;<\/span>
\nf<\/span>(<\/span>t +<\/span> c5 *<\/span> h,<\/span> tmp,<\/span> k5)<\/span>;<\/span><\/p>\n

\/\/ k6<\/span>
\nfor<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> n;<\/span> ++<\/span>i)<\/span>
\n tmp[<\/span>i]<\/span> =<\/span> y[<\/span>i]<\/span> +<\/span> h *<\/span> (<\/span>a61 *<\/span> k1[<\/span>i]<\/span> +<\/span> a62 *<\/span> k2[<\/span>i]<\/span> +<\/span> a63 *<\/span> k3[<\/span>i]<\/span> +<\/span> a64 *<\/span> k4[<\/span>i]<\/span> +<\/span> a65 *<\/span> k5[<\/span>i]<\/span>)<\/span>;<\/span>
\nf<\/span>(<\/span>t +<\/span> c6 *<\/span> h,<\/span> tmp,<\/span> k6)<\/span>;<\/span><\/p>\n

\/\/ 5th order solution y_next using b1..b6<\/span>
\ny_next.<\/span>resize<\/span>(<\/span>n)<\/span>;<\/span>
\nfor<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> n;<\/span> ++<\/span>i)<\/span>
\n y_next[<\/span>i]<\/span> =<\/span> y[<\/span>i]<\/span> +<\/span> h *<\/span> (<\/span>b1 *<\/span> k1[<\/span>i]<\/span> +<\/span> b2 *<\/span> k2[<\/span>i]<\/span> +<\/span> b3 *<\/span> k3[<\/span>i]<\/span> +<\/span> b4 *<\/span> k4[<\/span>i]<\/span> +<\/span> b5 *<\/span> k5[<\/span>i]<\/span> +<\/span> b6 *<\/span> k6[<\/span>i]<\/span> +<\/span> b7 *<\/span> 0.0<\/span>)<\/span>;<\/span><\/p>\n

\/\/ k7 evaluated at (t+h, y_next)<\/span>
\nf<\/span>(<\/span>t +<\/span> c7 *<\/span> h,<\/span> y_next,<\/span> k7)<\/span>;<\/span><\/p>\n

\/\/ error estimate (difference between 5th and 4th order solutions)<\/span>
\nerr.<\/span>resize<\/span>(<\/span>n)<\/span>;<\/span>
\nfor<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> n;<\/span> ++<\/span>i)<\/span> {<\/span>
\n const<\/span> double<\/span> e =<\/span> h *<\/span> (<\/span>(<\/span>b1 –<\/span> bh1)<\/span> *<\/span> k1[<\/span>i]<\/span> +<\/span> (<\/span>b2 –<\/span> bh2)<\/span> *<\/span> k2[<\/span>i]<\/span> +<\/span> (<\/span>b3 –<\/span> bh3)<\/span> *<\/span> k3[<\/span>i]<\/span> +<\/span> (<\/span>b4 –<\/span> bh4)<\/span> *<\/span> k4[<\/span>i]<\/span> +<\/span>
\n (<\/span>b5 –<\/span> bh5)<\/span> *<\/span> k5[<\/span>i]<\/span> +<\/span> (<\/span>b6 –<\/span> bh6)<\/span> *<\/span> k6[<\/span>i]<\/span> +<\/span> (<\/span>0.0<\/span> –<\/span> bh7)<\/span> *<\/span> k7[<\/span>i]<\/span>)<\/span>;<\/span>
\n err[<\/span>i]<\/span> =<\/span> e;<\/span>
\n}<\/span><\/p>\n

\/\/ Export stages for dense output interpolation<\/span>
\nK[<\/span>0<\/span>]<\/span> =<\/span> std::<\/span>move<\/span>(<\/span>k1)<\/span>;<\/span>
\nK[<\/span>1<\/span>]<\/span> =<\/span> std::<\/span>move<\/span>(<\/span>k2)<\/span>;<\/span>
\nK[<\/span>2<\/span>]<\/span> =<\/span> std::<\/span>move<\/span>(<\/span>k3)<\/span>;<\/span>
\nK[<\/span>3<\/span>]<\/span> =<\/span> std::<\/span>move<\/span>(<\/span>k4)<\/span>;<\/span>
\nK[<\/span>4<\/span>]<\/span> =<\/span> std::<\/span>move<\/span>(<\/span>k5)<\/span>;<\/span>
\nK[<\/span>5<\/span>]<\/span> =<\/span> std::<\/span>move<\/span>(<\/span>k6)<\/span>;<\/span>
\nK[<\/span>6<\/span>]<\/span> =<\/span> std::<\/span>move<\/span>(<\/span>k7)<\/span>;<\/span>
\n}<\/span><\/p>\n

\/\/ Adaptive integration wrapper.<\/span>
\nstatic<\/span> State integrate_runge_kutta_dopri5<\/span>(<\/span>const<\/span> DerivFunc&<\/span> f,<\/span>
\n double<\/span> t0,<\/span>
\n const<\/span> State&<\/span> y0,<\/span>
\n double<\/span> t1,<\/span>
\n double<\/span> h_init,<\/span>
\n double<\/span> atol,<\/span>
\n double<\/span> rtol,<\/span>
\n bool<\/span> normControl =<\/span> false<\/span>,<\/span>
\n int<\/span> max_steps =<\/span> 100000<\/span>)<\/span> {<\/span>
\n if<\/span> (<\/span>h_init ==<\/span> 0.0<\/span>)<\/span> {<\/span>
\n throw<\/span> std::<\/span>runtime_error<\/span>(<\/span>"h_init must be non-zero"<\/span>)<\/span>;<\/span>
\n }<\/span>
\n if<\/span> (<\/span>max_steps <=<\/span> 0<\/span>)<\/span> {<\/span>
\n throw<\/span> std::<\/span>runtime_error<\/span>(<\/span>"max_steps must be positive"<\/span>)<\/span>;<\/span>
\n }<\/span><\/p>\n

const<\/span> double<\/span> dir =<\/span> (<\/span>t1 >=<\/span> t0)<\/span> ?<\/span> 1.0<\/span> :<\/span> –<\/span>1.0<\/span>;<\/span>
\n double<\/span> t =<\/span> t0;<\/span>
\n State y =<\/span> y0;<\/span>
\n double<\/span> h =<\/span> std::<\/span>abs<\/span>(<\/span>h_init)<\/span> *<\/span> dir;<\/span><\/p>\n

\/\/ Step-size control parameters<\/span>
\n const<\/span> double<\/span> safety =<\/span> 0.9<\/span>;<\/span>
\n const<\/span> double<\/span> min_factor =<\/span> 0.2<\/span>;<\/span>
\n const<\/span> double<\/span> max_factor =<\/span> 10.0<\/span>;<\/span>
\n const<\/span> double<\/span> tiny =<\/span> std::<\/span>numeric_limits<\/span><<\/span>double<\/span>><\/span>::<\/span>min<\/span>(<\/span>)<\/span>;<\/span><\/p>\n

for<\/span> (<\/span>int<\/span> step =<\/span> 0<\/span>;<\/span> step <<\/span> max_steps;<\/span> ++<\/span>step)<\/span> {<\/span>
\n if<\/span> (<\/span>(<\/span>dir ><\/span> 0.0<\/span> &&<\/span> t >=<\/span> t1)<\/span> ||<\/span> (<\/span>dir <<\/span> 0.0<\/span> &&<\/span> t <=<\/span> t1)<\/span>)<\/span> {<\/span>
\n return<\/span> y;<\/span>
\n }<\/span>
\n \/\/ Clamp step to hit t1 exactly<\/span>
\n if<\/span> (<\/span>(<\/span>dir ><\/span> 0.0<\/span> &&<\/span> t +<\/span> h ><\/span> t1)<\/span> ||<\/span> (<\/span>dir <<\/span> 0.0<\/span> &&<\/span> t +<\/span> h <<\/span> t1)<\/span>)<\/span> {<\/span>
\n h =<\/span> t1 –<\/span> t;<\/span>
\n }<\/span><\/p>\n

State y_next,<\/span> err;<\/span>
\n std::<\/span>array<<\/span>State,<\/span>7<\/span>><\/span> K;<\/span>
\n runge_kutta_dopri5_step<\/span>(<\/span>f,<\/span> t,<\/span> y,<\/span> h,<\/span> y_next,<\/span> err,<\/span> K)<\/span>;<\/span>
\n const<\/span> double<\/span> err_norm =<\/span> normControl
\n ?<\/span> two_norm_scaled<\/span>(<\/span>err,<\/span> y,<\/span> y_next,<\/span> atol,<\/span> rtol)<\/span>
\n :<\/span> inf_norm_scaled<\/span>(<\/span>err,<\/span> y,<\/span> y_next,<\/span> atol,<\/span> rtol)<\/span>;<\/span><\/p>\n

if<\/span> (<\/span>err_norm <=<\/span> 1.0<\/span>)<\/span> {<\/span>
\n \/\/ Accept step<\/span>
\n t +=<\/span> h;<\/span>
\n y.<\/span>swap<\/span>(<\/span>y_next)<\/span>;<\/span>
\n }<\/span><\/p>\n

\/\/ Compute next step size<\/span>
\n \/\/ Embedded RK 5(4): controller exponent 1\/5.<\/span>
\n double<\/span> factor;<\/span>
\n if<\/span> (<\/span>err_norm <=<\/span> tiny)<\/span> {<\/span>
\n factor =<\/span> max_factor;<\/span>
\n }<\/span> else<\/span> if<\/span> (<\/span>!<\/span>std::<\/span>isfinite<\/span>(<\/span>err_norm)<\/span>)<\/span> {<\/span>
\n \/\/ Non-finite error: shrink step aggressively<\/span>
\n factor =<\/span> min_factor;<\/span>
\n }<\/span> else<\/span> {<\/span>
\n factor =<\/span> safety *<\/span> std::<\/span>pow<\/span>(<\/span>1.0<\/span> \/<\/span> err_norm,<\/span> 0.2<\/span>)<\/span>;<\/span><\/p>\n

factor =<\/span> std::<\/span>min<\/span>(<\/span>max_factor,<\/span> std::<\/span>max<\/span>(<\/span>min_factor,<\/span> factor)<\/span>)<\/span>;<\/span>
\n}<\/span>
\n h *=<\/span> factor;<\/span>
\n if<\/span> (<\/span>h ==<\/span> 0.0<\/span>)<\/span> {<\/span>
\n throw<\/span> std::<\/span>runtime_error<\/span>(<\/span>"Step size underflow"<\/span>)<\/span>;<\/span>
\n }<\/span>
\n \/\/ If rejected (including non-finite), try again without advancing time\/state.<\/span>
\n if<\/span> (<\/span>err_norm ><\/span> 1.0<\/span> ||<\/span> !<\/span>std::<\/span>isfinite<\/span>(<\/span>err_norm)<\/span>)<\/span> {<\/span>
\n continue<\/span>;<\/span>
\n }<\/span>
\n }<\/span><\/p>\n

throw<\/span> std::<\/span>runtime_error<\/span>(<\/span>"Reached max_steps without finishing integration"<\/span>)<\/span>;<\/span>
\n}<\/span><\/p>\n

struct<\/span> Trajectory<\/span> {<\/span>
\n std::<\/span>vector<<\/span>double<\/span>><\/span> t;<\/span>
\n std::<\/span>vector<<\/span>State><\/span> y;<\/span>
\n}<\/span>;<\/span><\/p>\n

\/\/ Integrate and return dense output samples according to refine (>=1)<\/span>
\nstatic<\/span> Trajectory integrate_runge_kutta_dopri5_refined<\/span>(<\/span>const<\/span> DerivFunc&<\/span> f,<\/span>
\n double<\/span> t0,<\/span>
\n const<\/span> State&<\/span> y0,<\/span>
\n double<\/span> t1,<\/span>
\n double<\/span> h_init,<\/span>
\n double<\/span> atol,<\/span>
\n double<\/span> rtol,<\/span>
\n bool<\/span> normControl,<\/span>
\n int<\/span> refine,<\/span>
\n int<\/span> max_steps =<\/span> 100000<\/span>)<\/span> {<\/span>
\n if<\/span> (<\/span>h_init ==<\/span> 0.0<\/span>)<\/span> throw<\/span> std::<\/span>runtime_error<\/span>(<\/span>"h_init must be non-zero"<\/span>)<\/span>;<\/span>
\n if<\/span> (<\/span>max_steps <=<\/span> 0<\/span>)<\/span> throw<\/span> std::<\/span>runtime_error<\/span>(<\/span>"max_steps must be positive"<\/span>)<\/span>;<\/span>
\n if<\/span> (<\/span>refine <<\/span> 1<\/span>)<\/span> refine =<\/span> 1<\/span>;<\/span><\/p>\n

const<\/span> double<\/span> dir =<\/span> (<\/span>t1 >=<\/span> t0)<\/span> ?<\/span> 1.0<\/span> :<\/span> –<\/span>1.0<\/span>;<\/span>
\n double<\/span> t =<\/span> t0;<\/span>
\n State y =<\/span> y0;<\/span>
\n double<\/span> h =<\/span> std::<\/span>abs<\/span>(<\/span>h_init)<\/span> *<\/span> dir;<\/span><\/p>\n

const<\/span> double<\/span> safety =<\/span> 0.9<\/span>;<\/span>
\n const<\/span> double<\/span> min_factor =<\/span> 0.2<\/span>;<\/span>
\n const<\/span> double<\/span> max_factor =<\/span> 10.0<\/span>;<\/span>
\n const<\/span> double<\/span> tiny =<\/span> std::<\/span>numeric_limits<\/span><<\/span>double<\/span>><\/span>::<\/span>min<\/span>(<\/span>)<\/span>;<\/span><\/p>\n

Trajectory traj;<\/span>
\n traj.<\/span>t.<\/span>push_back<\/span>(<\/span>t)<\/span>;<\/span>
\n traj.<\/span>y.<\/span>push_back<\/span>(<\/span>y)<\/span>;<\/span><\/p>\n

for<\/span> (<\/span>int<\/span> step =<\/span> 0<\/span>;<\/span> step <<\/span> max_steps;<\/span> ++<\/span>step)<\/span> {<\/span>
\n if<\/span> (<\/span>(<\/span>dir ><\/span> 0.0<\/span> &&<\/span> t >=<\/span> t1)<\/span> ||<\/span> (<\/span>dir <<\/span> 0.0<\/span> &&<\/span> t <=<\/span> t1)<\/span>)<\/span> {<\/span>
\n return<\/span> traj;<\/span>
\n }<\/span><\/p>\n

if<\/span> (<\/span>(<\/span>dir ><\/span> 0.0<\/span> &&<\/span> t +<\/span> h ><\/span> t1)<\/span> ||<\/span> (<\/span>dir <<\/span> 0.0<\/span> &&<\/span> t +<\/span> h <<\/span> t1)<\/span>)<\/span> {<\/span>
\n h =<\/span> t1 –<\/span> t;<\/span>
\n }<\/span><\/p>\n

State y_next,<\/span> err;<\/span>
\n std::<\/span>array<<\/span>State,<\/span>7<\/span>><\/span> K;<\/span>
\n runge_kutta_dopri5_step<\/span>(<\/span>f,<\/span> t,<\/span> y,<\/span> h,<\/span> y_next,<\/span> err,<\/span> K)<\/span>;<\/span>
\n const<\/span> double<\/span> err_norm =<\/span> normControl
\n ?<\/span> two_norm_scaled<\/span>(<\/span>err,<\/span> y,<\/span> y_next,<\/span> atol,<\/span> rtol)<\/span>
\n :<\/span> inf_norm_scaled<\/span>(<\/span>err,<\/span> y,<\/span> y_next,<\/span> atol,<\/span> rtol)<\/span>;<\/span><\/p>\n

if<\/span> (<\/span>err_norm <=<\/span> 1.0<\/span>)<\/span> {<\/span>
\n \/\/ Accepted: emit interior refined points using dense output<\/span>
\n const<\/span> double<\/span> dt =<\/span> h;<\/span>
\n for<\/span> (<\/span>int<\/span> i =<\/span> 1<\/span>;<\/span> i <<\/span> refine;<\/span> ++<\/span>i)<\/span> {<\/span>
\n const<\/span> double<\/span> s =<\/span> static_cast<\/span><<\/span>double<\/span>><\/span><\/span><\/span>(<\/span>i)<\/span> \/<\/span> static_cast<\/span><<\/span>double<\/span>><\/span><\/span><\/span>(<\/span>refine)<\/span>;<\/span>
\n State yi =<\/span> dopri5_dense_output<\/span>(<\/span>dt,<\/span> y,<\/span> y_next,<\/span> K,<\/span> s)<\/span>;<\/span>
\n traj.<\/span>t.<\/span>push_back<\/span>(<\/span>t +<\/span> s *<\/span> h)<\/span>;<\/span>
\n traj.<\/span>y.<\/span>push_back<\/span>(<\/span>std::<\/span>move<\/span>(<\/span>yi)<\/span>)<\/span>;<\/span>
\n }<\/span>
\n \/\/ Endpoint<\/span>
\n t +=<\/span> h;<\/span>
\n y.<\/span>swap<\/span>(<\/span>y_next)<\/span>;<\/span>
\n traj.<\/span>t.<\/span>push_back<\/span>(<\/span>t)<\/span>;<\/span>
\n traj.<\/span>y.<\/span>push_back<\/span>(<\/span>y)<\/span>;<\/span>
\n }<\/span><\/p>\n

\/\/ Step-size controller<\/span>
\n double<\/span> factor;<\/span>
\n if<\/span> (<\/span>err_norm <=<\/span> tiny)<\/span> {<\/span>
\n factor =<\/span> max_factor;<\/span>
\n }<\/span> else<\/span> if<\/span> (<\/span>!<\/span>std::<\/span>isfinite<\/span>(<\/span>err_norm)<\/span>)<\/span> {<\/span>
\n factor =<\/span> min_factor;<\/span>
\n }<\/span> else<\/span> {<\/span>
\n factor =<\/span> safety *<\/span> std::<\/span>pow<\/span>(<\/span>1.0<\/span> \/<\/span> err_norm,<\/span> 0.2<\/span>)<\/span>;<\/span>
\n factor =<\/span> std::<\/span>min<\/span>(<\/span>max_factor,<\/span> std::<\/span>max<\/span>(<\/span>min_factor,<\/span> factor)<\/span>)<\/span>;<\/span>
\n }<\/span><\/p>\n

h *=<\/span> factor;<\/span>
\n if<\/span> (<\/span>h ==<\/span> 0.0<\/span>)<\/span> throw<\/span> std::<\/span>runtime_error<\/span>(<\/span>"Step size underflow"<\/span>)<\/span>;<\/span>
\n if<\/span> (<\/span>err_norm ><\/span> 1.0<\/span> ||<\/span> !<\/span>std::<\/span>isfinite<\/span>(<\/span>err_norm)<\/span>)<\/span> continue<\/span>;<\/span>
\n }<\/span><\/p>\n

throw<\/span> std::<\/span>runtime_error<\/span>(<\/span>"Reached max_steps without finishing integration"<\/span>)<\/span>;<\/span>
\n}<\/span><\/p>\n

\/\/ Integrate and return samples exactly at requested times using dense output per accepted step.<\/span>
\nstatic<\/span> Trajectory integrate_runge_kutta_dopri5_at_times<\/span>(<\/span>const<\/span> DerivFunc&<\/span> f,<\/span>
\n double<\/span> t0,<\/span>
\n const<\/span> State&<\/span> y0,<\/span>
\n const<\/span> std::<\/span>vector<<\/span>double<\/span>><\/span>&<\/span> times,<\/span>
\n double<\/span> h_init,<\/span>
\n double<\/span> atol,<\/span>
\n double<\/span> rtol,<\/span>
\n bool<\/span> normControl,<\/span>
\n int<\/span> max_steps =<\/span> 100000<\/span>)<\/span> {<\/span>
\n if<\/span> (<\/span>h_init ==<\/span> 0.0<\/span>)<\/span> throw<\/span> std::<\/span>runtime_error<\/span>(<\/span>"h_init must be non-zero"<\/span>)<\/span>;<\/span>
\n if<\/span> (<\/span>max_steps <=<\/span> 0<\/span>)<\/span> throw<\/span> std::<\/span>runtime_error<\/span>(<\/span>"max_steps must be positive"<\/span>)<\/span>;<\/span>
\n if<\/span> (<\/span>times.<\/span>empty<\/span>(<\/span>)<\/span>)<\/span> throw<\/span> std::<\/span>runtime_error<\/span>(<\/span>"times must be non-empty"<\/span>)<\/span>;<\/span>
\n \/\/ Ensure monotonicity direction matches integration<\/span>
\n const<\/span> double<\/span> dir =<\/span> (<\/span>times.<\/span>back<\/span>(<\/span>)<\/span> >=<\/span> times.<\/span>front<\/span>(<\/span>)<\/span>)<\/span> ?<\/span> 1.0<\/span> :<\/span> –<\/span>1.0<\/span>;<\/span>
\n \/\/ t0 must equal first time to match initial condition semantics<\/span>
\n if<\/span> (<\/span>std::<\/span>abs<\/span>(<\/span>times.<\/span>front<\/span>(<\/span>)<\/span> –<\/span> t0)<\/span> ><\/span> 0<\/span>)<\/span> {<\/span>
\n t0 =<\/span> times.<\/span>front<\/span>(<\/span>)<\/span>;<\/span>
\n }<\/span><\/p>\n

double<\/span> t =<\/span> t0;<\/span>
\n State y =<\/span> y0;<\/span>
\n double<\/span> h =<\/span> std::<\/span>abs<\/span>(<\/span>h_init)<\/span> *<\/span> dir;<\/span><\/p>\n

const<\/span> double<\/span> safety =<\/span> 0.9<\/span>;<\/span>
\n const<\/span> double<\/span> min_factor =<\/span> 0.2<\/span>;<\/span>
\n const<\/span> double<\/span> max_factor =<\/span> 10.0<\/span>;<\/span>
\n const<\/span> double<\/span> tiny =<\/span> std::<\/span>numeric_limits<\/span><<\/span>double<\/span>><\/span>::<\/span>min<\/span>(<\/span>)<\/span>;<\/span><\/p>\n

Trajectory traj;<\/span>
\n size_t idx =<\/span> 0<\/span>;<\/span>
\n \/\/ Emit initial sample if requested equals t0<\/span>
\n if<\/span> (<\/span>(<\/span>dir ><\/span> 0<\/span> &&<\/span> times[<\/span>idx]<\/span> <=<\/span> t)<\/span> ||<\/span> (<\/span>dir <<\/span> 0<\/span> &&<\/span> times[<\/span>idx]<\/span> >=<\/span> t)<\/span>)<\/span> {<\/span>
\n traj.<\/span>t.<\/span>push_back<\/span>(<\/span>times[<\/span>idx]<\/span>)<\/span>;<\/span>
\n traj.<\/span>y.<\/span>push_back<\/span>(<\/span>y)<\/span>;<\/span>
\n ++<\/span>idx;<\/span>
\n }<\/span><\/p>\n

for<\/span> (<\/span>int<\/span> step =<\/span> 0<\/span>;<\/span> step <<\/span> max_steps &&<\/span> idx <<\/span> times.<\/span>size<\/span>(<\/span>)<\/span>;<\/span> ++<\/span>step)<\/span> {<\/span>
\n \/\/ If we've passed the last requested time, stop<\/span>
\n if<\/span> (<\/span>(<\/span>dir ><\/span> 0.0<\/span> &&<\/span> t >=<\/span> times.<\/span>back<\/span>(<\/span>)<\/span>)<\/span> ||<\/span> (<\/span>dir <<\/span> 0.0<\/span> &&<\/span> t <=<\/span> times.<\/span>back<\/span>(<\/span>)<\/span>)<\/span>)<\/span> {<\/span>
\n break<\/span>;<\/span>
\n }<\/span><\/p>\n

\/\/ Adjust step to not skip beyond last requested time<\/span>
\n if<\/span> (<\/span>(<\/span>dir ><\/span> 0.0<\/span> &&<\/span> t +<\/span> h ><\/span> times.<\/span>back<\/span>(<\/span>)<\/span>)<\/span> ||<\/span> (<\/span>dir <<\/span> 0.0<\/span> &&<\/span> t +<\/span> h <<\/span> times.<\/span>back<\/span>(<\/span>)<\/span>)<\/span>)<\/span> {<\/span>
\n h =<\/span> times.<\/span>back<\/span>(<\/span>)<\/span> –<\/span> t;<\/span>
\n }<\/span><\/p>\n

State y_next,<\/span> err;<\/span>
\n std::<\/span>array<<\/span>State,<\/span>7<\/span>><\/span> K;<\/span>
\n runge_kutta_dopri5_step<\/span>(<\/span>f,<\/span> t,<\/span> y,<\/span> h,<\/span> y_next,<\/span> err,<\/span> K)<\/span>;<\/span>
\n const<\/span> double<\/span> err_norm =<\/span> normControl
\n ?<\/span> two_norm_scaled<\/span>(<\/span>err,<\/span> y,<\/span> y_next,<\/span> atol,<\/span> rtol)<\/span>
\n :<\/span> inf_norm_scaled<\/span>(<\/span>err,<\/span> y,<\/span> y_next,<\/span> atol,<\/span> rtol)<\/span>;<\/span><\/p>\n

\/\/ Step-size controller<\/span>
\n double<\/span> factor;<\/span>
\n if<\/span> (<\/span>err_norm <=<\/span> tiny)<\/span> {<\/span>
\n factor =<\/span> max_factor;<\/span>
\n }<\/span> else<\/span> if<\/span> (<\/span>!<\/span>std::<\/span>isfinite<\/span>(<\/span>err_norm)<\/span>)<\/span> {<\/span>
\n factor =<\/span> min_factor;<\/span>
\n }<\/span> else<\/span> {<\/span>
\n factor =<\/span> safety *<\/span> std::<\/span>pow<\/span>(<\/span>1.0<\/span> \/<\/span> err_norm,<\/span> 0.2<\/span>)<\/span>;<\/span>
\n factor =<\/span> std::<\/span>min<\/span>(<\/span>max_factor,<\/span> std::<\/span>max<\/span>(<\/span>min_factor,<\/span> factor)<\/span>)<\/span>;<\/span>
\n }<\/span>
\n h *=<\/span> factor;<\/span>
\n if<\/span> (<\/span>h ==<\/span> 0.0<\/span>)<\/span> throw<\/span> std::<\/span>runtime_error<\/span>(<\/span>"Step size underflow"<\/span>)<\/span>;<\/span>
\n if<\/span> (<\/span>err_norm ><\/span> 1.0<\/span> ||<\/span> !<\/span>std::<\/span>isfinite<\/span>(<\/span>err_norm)<\/span>)<\/span> continue<\/span>;<\/span> \/\/ rejected<\/span><\/p>\n

\/\/ Accepted: emit samples at times within (t, t+h]<\/span>
\n const<\/span> double<\/span> dt =<\/span> times.<\/span>back<\/span>(<\/span>)<\/span> –<\/span> times.<\/span>front<\/span>(<\/span>)<\/span>;<\/span> \/\/ unused; keep h for s calc<\/span>
\n while<\/span> (<\/span>idx <<\/span> times.<\/span>size<\/span>(<\/span>)<\/span>)<\/span> {<\/span>
\n const<\/span> double<\/span> tau =<\/span> times[<\/span>idx]<\/span>;<\/span>
\n \/\/ Check whether tau is within this step interval<\/span>
\n const<\/span> bool<\/span> in_interval =<\/span> (<\/span>dir ><\/span> 0.0<\/span>)<\/span> ?<\/span> (<\/span>tau ><\/span> t &&<\/span> tau <=<\/span> t +<\/span> h)<\/span> :<\/span> (<\/span>tau <<\/span> t &&<\/span> tau >=<\/span> t +<\/span> h)<\/span>;<\/span>
\n if<\/span> (<\/span>!<\/span>in_interval)<\/span> break<\/span>;<\/span>
\n const<\/span> double<\/span> s =<\/span> (<\/span>tau –<\/span> t)<\/span> \/<\/span> h;<\/span> \/\/ works both directions<\/span>
\n State yi =<\/span> dopri5_dense_output<\/span>(<\/span>std::<\/span>abs<\/span>(<\/span>h)<\/span>,<\/span> y,<\/span> y_next,<\/span> K,<\/span> s)<\/span>;<\/span>
\n traj.<\/span>t.<\/span>push_back<\/span>(<\/span>tau)<\/span>;<\/span>
\n traj.<\/span>y.<\/span>push_back<\/span>(<\/span>std::<\/span>move<\/span>(<\/span>yi)<\/span>)<\/span>;<\/span>
\n ++<\/span>idx;<\/span>
\n }<\/span><\/p>\n

\/\/ Advance to end of step<\/span>
\n t +=<\/span> h;<\/span>
\n y.<\/span>swap<\/span>(<\/span>y_next)<\/span>;<\/span>
\n }<\/span><\/p>\n

if<\/span> (<\/span>idx !=<\/span> times.<\/span>size<\/span>(<\/span>)<\/span>)<\/span> {<\/span>
\n throw<\/span> std::<\/span>runtime_error<\/span>(<\/span>"Failed to sample all requested times before reaching max_steps"<\/span>)<\/span>;<\/span>
\n }<\/span>
\n return<\/span> traj;<\/span>
\n}<\/span><\/p>\n

struct<\/span> Options<\/span> {<\/span>
\n double<\/span> t0 =<\/span> 0.0<\/span>;<\/span>
\n double<\/span> t1 =<\/span> 5.0<\/span>;<\/span>
\n double<\/span> h_init =<\/span> 0.1<\/span>;<\/span>
\n double<\/span> atol =<\/span> 1e-10<\/span>;<\/span>
\n double<\/span> rtol =<\/span> 1e-10<\/span>;<\/span>
\n bool<\/span> normControl =<\/span> false<\/span>;<\/span> \/\/ Octave default is off<\/span>
\n int<\/span> refine =<\/span> 4<\/span>;<\/span> \/\/ Octave default is 4 for ode45<\/span>
\n int<\/span> max_steps =<\/span> 100000<\/span>;<\/span>
\n std::<\/span>vector<<\/span>double<\/span>><\/span> times;<\/span> \/\/ explicit sample times; overrides refine when non-empty<\/span>
\n}<\/span>;<\/span><\/p>\n

static<\/span> Options parseArgs<\/span>(<\/span>int<\/span> argc,<\/span> char<\/span>*<\/span>*<\/span> argv,<\/span> const<\/span> Options&<\/span> defaults)<\/span> {<\/span>
\n Options opt =<\/span> defaults;<\/span>
\n auto<\/span> starts_with =<\/span> [<\/span>]<\/span>(<\/span>const<\/span> std::<\/span>string&<\/span> s,<\/span> const<\/span> std::<\/span>string&<\/span> p)<\/span> {<\/span> return<\/span> s.<\/span>rfind<\/span>(<\/span>p,<\/span>0<\/span>)<\/span>==<\/span>0<\/span>;<\/span> }<\/span>;<\/span>
\n for<\/span> (<\/span>int<\/span> i =<\/span> 1<\/span>;<\/span> i <<\/span> argc;<\/span> ++<\/span>i)<\/span> {<\/span>
\n std::<\/span>string a =<\/span> argv[<\/span>i]<\/span>;<\/span>
\n if<\/span> (<\/span>a ==<\/span> "–help"<\/span> ||<\/span> a ==<\/span> "-h"<\/span>)<\/span> {<\/span>
\n std::<\/span>cout <<<\/span> "Options:\\\\n"<\/span>
\n <<<\/span> " –t0=VAL start time (default 0)\\\\n"<\/span>
\n <<<\/span> " –t1=VAL end time (default 5)\\\\n"<\/span>
\n <<<\/span> " –h-init=VAL initial step size\\\\n"<\/span>
\n <<<\/span> " –atol=VAL absolute tolerance\\\\n"<\/span>
\n <<<\/span> " –rtol=VAL relative tolerance\\\\n"<\/span>
\n <<<\/span> " –normcontrol=on|off norm control (default off)\\\\n"<\/span>
\n <<<\/span> " –refine=N points per step (default 4)\\\\n"<\/span>
\n <<<\/span> " –maxsteps=N max steps (default 100000)\\\\n"<\/span>
\n <<<\/span> " –times=t0,t1,… explicit times to sample (overrides refine)\\\\n"<\/span>;<\/span>
\n std::<\/span>exit<\/span>(<\/span>0<\/span>)<\/span>;<\/span>
\n }<\/span>
\n auto<\/span> eq =<\/span> a.<\/span>find<\/span>(<\/span>'='<\/span>)<\/span>;<\/span>
\n if<\/span> (<\/span>eq ==<\/span> std::<\/span>string::<\/span>npos)<\/span> continue<\/span>;<\/span>
\n const<\/span> std::<\/span>string key =<\/span> a.<\/span>substr<\/span>(<\/span>0<\/span>,<\/span> eq+<\/span>1<\/span>)<\/span>;<\/span>
\n const<\/span> std::<\/span>string val =<\/span> a.<\/span>substr<\/span>(<\/span>eq+<\/span>1<\/span>)<\/span>;<\/span>
\n try<\/span> {<\/span>
\n if<\/span> (<\/span>key ==<\/span> "–t0="<\/span>)<\/span> opt.<\/span>t0 =<\/span> std::<\/span>stod<\/span>(<\/span>val)<\/span>;<\/span>
\n else<\/span> if<\/span> (<\/span>key ==<\/span> "–t1="<\/span>)<\/span> opt.<\/span>t1 =<\/span> std::<\/span>stod<\/span>(<\/span>val)<\/span>;<\/span>
\n else<\/span> if<\/span> (<\/span>key ==<\/span> "–h-init="<\/span>)<\/span> opt.<\/span>h_init =<\/span> std::<\/span>stod<\/span>(<\/span>val)<\/span>;<\/span>
\n else<\/span> if<\/span> (<\/span>key ==<\/span> "–atol="<\/span>)<\/span> opt.<\/span>atol =<\/span> std::<\/span>stod<\/span>(<\/span>val)<\/span>;<\/span>
\n else<\/span> if<\/span> (<\/span>key ==<\/span> "–rtol="<\/span>)<\/span> opt.<\/span>rtol =<\/span> std::<\/span>stod<\/span>(<\/span>val)<\/span>;<\/span>
\n else<\/span> if<\/span> (<\/span>key ==<\/span> "–refine="<\/span>)<\/span> opt.<\/span>refine =<\/span> std::<\/span>stoi<\/span>(<\/span>val)<\/span>;<\/span>
\n else<\/span> if<\/span> (<\/span>key ==<\/span> "–maxsteps="<\/span>)<\/span> opt.<\/span>max_steps =<\/span> std::<\/span>stoi<\/span>(<\/span>val)<\/span>;<\/span>
\n else<\/span> if<\/span> (<\/span>key ==<\/span> "–normcontrol="<\/span>)<\/span> {<\/span>
\n const<\/span> std::<\/span>string v =<\/span> val;<\/span>
\n opt.<\/span>normControl =<\/span> (<\/span>v ==<\/span> "on"<\/span> ||<\/span> v ==<\/span> "ON"<\/span> ||<\/span> v ==<\/span> "true"<\/span> ||<\/span> v ==<\/span> "1"<\/span>)<\/span>;<\/span>
\n }<\/span> else<\/span> if<\/span> (<\/span>key ==<\/span> "–times="<\/span>)<\/span> {<\/span>
\n \/\/ Parse comma-separated list of doubles<\/span>
\n opt.<\/span>times.<\/span>clear<\/span>(<\/span>)<\/span>;<\/span>
\n size_t start =<\/span> 0<\/span>;<\/span> while<\/span> (<\/span>start <=<\/span> val.<\/span>size<\/span>(<\/span>)<\/span>)<\/span> {<\/span>
\n size_t comma =<\/span> val.<\/span>find<\/span>(<\/span>','<\/span>,<\/span> start)<\/span>;<\/span>
\n std::<\/span>string tok =<\/span> val.<\/span>substr<\/span>(<\/span>start,<\/span> (<\/span>comma==<\/span>std::<\/span>string::<\/span>npos)<\/span>?<\/span> std::<\/span>string::<\/span>npos :<\/span> (<\/span>comma–<\/span>start)<\/span>)<\/span>;<\/span>
\n if<\/span> (<\/span>!<\/span>tok.<\/span>empty<\/span>(<\/span>)<\/span>)<\/span> opt.<\/span>times.<\/span>push_back<\/span>(<\/span>std::<\/span>stod<\/span>(<\/span>tok)<\/span>)<\/span>;<\/span>
\n if<\/span> (<\/span>comma ==<\/span> std::<\/span>string::<\/span>npos)<\/span> break<\/span>;<\/span> else<\/span> start =<\/span> comma +<\/span> 1<\/span>;<\/span>
\n }<\/span>
\n }<\/span>
\n }<\/span> catch<\/span> (<\/span>.<\/span>.<\/span>.<\/span>)<\/span> {<\/span>
\n std::<\/span>cerr <<<\/span> "Invalid value for option: "<\/span> <<<\/span> a <<<\/span> "\\\\n"<\/span>;<\/span>
\n std::<\/span>exit<\/span>(<\/span>1<\/span>)<\/span>;<\/span>
\n }<\/span>
\n }<\/span>
\n if<\/span> (<\/span>opt.<\/span>refine <<\/span> 1<\/span>)<\/span> opt.<\/span>refine =<\/span> 1<\/span>;<\/span>
\n return<\/span> opt;<\/span>
\n}<\/span><\/p>\n

int<\/span> main<\/span>(<\/span>int<\/span> argc,<\/span> char<\/span>*<\/span>*<\/span> argv)<\/span> {<\/span>
\n \/\/ Example: y' = -y, y(0) = 1. Exact solution: y(t) = exp(-t)<\/span>
\n DerivFunc f =<\/span> [<\/span>]<\/span>(<\/span>double<\/span> \/*t*\/<\/span>,<\/span> const<\/span> State&<\/span> y,<\/span> State&<\/span> dydt)<\/span> {<\/span>
\n dydt.<\/span>resize<\/span>(<\/span>y.<\/span>size<\/span>(<\/span>)<\/span>)<\/span>;<\/span>
\n for<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> y.<\/span>size<\/span>(<\/span>)<\/span>;<\/span> ++<\/span>i)<\/span> {<\/span>
\n dydt[<\/span>i]<\/span> =<\/span> –<\/span>y[<\/span>i]<\/span>;<\/span>
\n }<\/span>
\n }<\/span>;<\/span><\/p>\n

const<\/span> State y0 =<\/span> {<\/span>1.0<\/span>}<\/span>;<\/span>
\n Options defaults;<\/span>
\n defaults.<\/span>normControl =<\/span> true<\/span>;<\/span> \/\/ demo default on; Octave is off<\/span>
\n Options opt =<\/span> parseArgs<\/span>(<\/span>argc,<\/span> argv,<\/span> defaults)<\/span>;<\/span><\/p>\n

try<\/span> {<\/span>
\n \/\/ If explicit times provided, sample there; otherwise run default and optional refine<\/span>
\n if<\/span> (<\/span>!<\/span>opt.<\/span>times.<\/span>empty<\/span>(<\/span>)<\/span>)<\/span> {<\/span>
\n \/\/ Expect first time equals opt.t0; if not, we start at first time.<\/span>
\n auto<\/span> traj =<\/span> integrate_runge_kutta_dopri5_at_times<\/span>(<\/span>f,<\/span> opt.<\/span>t0,<\/span> y0,<\/span> opt.<\/span>times,<\/span> opt.<\/span>h_init,<\/span>
\n opt.<\/span>atol,<\/span> opt.<\/span>rtol,<\/span> opt.<\/span>normControl,<\/span> opt.<\/span>max_steps)<\/span>;<\/span>
\n std::<\/span>cout <<<\/span> "samples (first 5 of explicit times):\\\\n"<\/span>;<\/span>
\n for<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> std::<\/span>min<\/span><<\/span>size_t><\/span><\/span><\/span>(<\/span>traj.<\/span>t.<\/span>size<\/span>(<\/span>)<\/span>,<\/span>5<\/span>)<\/span>;<\/span> ++<\/span>i)<\/span> {<\/span>
\n std::<\/span>cout <<<\/span> std::<\/span>fixed <<<\/span> std::<\/span>setprecision<\/span>(<\/span>8<\/span>)<\/span> <<<\/span> traj.<\/span>t[<\/span>i]<\/span> <<<\/span> ": "<\/span> <<<\/span> traj.<\/span>y[<\/span>i]<\/span>[<\/span>0<\/span>]<\/span> <<<\/span> "\\\\n"<\/span>;<\/span>
\n }<\/span>
\n \/\/ Also print last sample<\/span>
\n std::<\/span>cout <<<\/span> std::<\/span>setprecision<\/span>(<\/span>16<\/span>)<\/span>;<\/span>
\n std::<\/span>cout <<<\/span> "y("<\/span> <<<\/span> traj.<\/span>t.<\/span>back<\/span>(<\/span>)<\/span> <<<\/span> ") = "<\/span> <<<\/span> traj.<\/span>y.<\/span>back<\/span>(<\/span>)<\/span>[<\/span>0<\/span>]<\/span> <<<\/span> "\\\\n"<\/span>;<\/span>
\n }<\/span> else<\/span> {<\/span>
\n const<\/span> State y1 =<\/span> integrate_runge_kutta_dopri5<\/span>(<\/span>f,<\/span> opt.<\/span>t0,<\/span> y0,<\/span> opt.<\/span>t1,<\/span> opt.<\/span>h_init,<\/span> opt.<\/span>atol,<\/span>
\n opt.<\/span>rtol,<\/span> opt.<\/span>normControl,<\/span> opt.<\/span>max_steps)<\/span>;<\/span>
\n const<\/span> double<\/span> exact =<\/span> std::<\/span>exp<\/span>(<\/span>–<\/span>opt.<\/span>t1)<\/span>;<\/span>
\n std::<\/span>cout <<<\/span> std::<\/span>setprecision<\/span>(<\/span>16<\/span>)<\/span>;<\/span>
\n std::<\/span>cout <<<\/span> "y("<\/span> <<<\/span> opt.<\/span>t1 <<<\/span> ") = "<\/span> <<<\/span> y1[<\/span>0<\/span>]<\/span> <<<\/span> "\\\\n"<\/span>;<\/span>
\n std::<\/span>cout <<<\/span> "exact = "<\/span> <<<\/span> exact <<<\/span> "\\\\n"<\/span>;<\/span>
\n std::<\/span>cout <<<\/span> "abs err = "<\/span> <<<\/span> std::<\/span>abs<\/span>(<\/span>y1[<\/span>0<\/span>]<\/span> –<\/span> exact)<\/span> <<<\/span> "\\\\n"<\/span>;<\/span><\/p>\n

if<\/span> (<\/span>opt.<\/span>refine ><\/span> 1<\/span>)<\/span> {<\/span>
\n auto<\/span> traj =<\/span> integrate_runge_kutta_dopri5_refined<\/span>(<\/span>f,<\/span> opt.<\/span>t0,<\/span> y0,<\/span> opt.<\/span>t1,<\/span> opt.<\/span>h_init,<\/span> opt.<\/span>atol,<\/span>
\n opt.<\/span>rtol,<\/span> opt.<\/span>normControl,<\/span> opt.<\/span>refine,<\/span> opt.<\/span>max_steps)<\/span>;<\/span>
\n std::<\/span>cout <<<\/span> "samples (first 5):\\\\n"<\/span>;<\/span>
\n for<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> std::<\/span>min<\/span><<\/span>size_t><\/span><\/span><\/span>(<\/span>traj.<\/span>t.<\/span>size<\/span>(<\/span>)<\/span>,<\/span>5<\/span>)<\/span>;<\/span> ++<\/span>i)<\/span> {<\/span>
\n std::<\/span>cout <<<\/span> std::<\/span>fixed <<<\/span> std::<\/span>setprecision<\/span>(<\/span>8<\/span>)<\/span> <<<\/span> traj.<\/span>t[<\/span>i]<\/span> <<<\/span> ": "<\/span> <<<\/span> traj.<\/span>y[<\/span>i]<\/span>[<\/span>0<\/span>]<\/span> <<<\/span> "\\\\n"<\/span>;<\/span>
\n }<\/span>
\n \/\/ Verify Octave's length relationship: len_R = R*len_1 – (R-1)<\/span>
\n auto<\/span> coarse =<\/span> integrate_runge_kutta_dopri5_refined<\/span>(<\/span>f,<\/span> opt.<\/span>t0,<\/span> y0,<\/span> opt.<\/span>t1,<\/span> opt.<\/span>h_init,<\/span> opt.<\/span>atol,<\/span>
\n opt.<\/span>rtol,<\/span> opt.<\/span>normControl,<\/span> 1<\/span>,<\/span> opt.<\/span>max_steps)<\/span>;<\/span>
\n const<\/span> size_t lenR =<\/span> traj.<\/span>t.<\/span>size<\/span>(<\/span>)<\/span>;<\/span>
\n const<\/span> size_t len1 =<\/span> coarse.<\/span>t.<\/span>size<\/span>(<\/span>)<\/span>;<\/span>
\n std::<\/span>cout <<<\/span> "length refine = "<\/span> <<<\/span> opt.<\/span>refine <<<\/span> ": "<\/span> <<<\/span> lenR <<<\/span> ", expected "<\/span> <<<\/span>
\n (<\/span>opt.<\/span>refine*<\/span>len1 –<\/span> (<\/span>opt.<\/span>refine–<\/span>1<\/span>)<\/span>)<\/span> <<<\/span> "\\\\n"<\/span>;<\/span>
\n }<\/span>
\n }<\/span>
\n return<\/span> 0<\/span>;<\/span>
\n }<\/span> catch<\/span> (<\/span>const<\/span> std::<\/span>exception&<\/span> e)<\/span> {<\/span>
\n std::<\/span>cerr <<<\/span> "Error: "<\/span> <<<\/span> e.<\/span>what<\/span>(<\/span>)<\/span> <<<\/span> "\\\\n"<\/span>;<\/span>
\n return<\/span> 1<\/span>;<\/span>
\n }<\/span>
\n}<\/span><\/p>\n

static<\/span> double<\/span> two_norm_scaled<\/span>(<\/span>const<\/span> State&<\/span> err,<\/span> const<\/span> State&<\/span> y,<\/span> const<\/span> State&<\/span> y_next,<\/span> double<\/span> atol,<\/span> double<\/span> rtol)<\/span> {<\/span>
\n \/\/ Matches two_norm_scaled with NormControl == true:<\/span>
\n \/\/ sc = max(AbsTol(:), RelTol * max(||y||_2, ||y_next||_2))<\/span>
\n \/\/ err_norm = ||err||_2 \/ sc<\/span>
\n if<\/span> (<\/span>err.<\/span>size<\/span>(<\/span>)<\/span> !=<\/span> y.<\/span>size<\/span>(<\/span>)<\/span> ||<\/span> err.<\/span>size<\/span>(<\/span>)<\/span> !=<\/span> y_next.<\/span>size<\/span>(<\/span>)<\/span>)<\/span> {<\/span>
\n throw<\/span> std::<\/span>runtime_error<\/span>(<\/span>"State size mismatch"<\/span>)<\/span>;<\/span>
\n }<\/span>
\n if<\/span> (<\/span>err.<\/span>empty<\/span>(<\/span>)<\/span>)<\/span> {<\/span>
\n return<\/span> 0.0<\/span>;<\/span>
\n }<\/span><\/p>\n

auto<\/span> norm2 =<\/span> [<\/span>]<\/span>(<\/span>const<\/span> State&<\/span> v)<\/span> {<\/span>
\n long<\/span> double<\/span> s =<\/span> 0.0L<\/span>;<\/span>
\n for<\/span> (<\/span>double<\/span> vi :<\/span> v)<\/span> s +=<\/span> static_cast<\/span><<\/span>long<\/span> double<\/span>><\/span><\/span><\/span>(<\/span>vi)<\/span> *<\/span> static_cast<\/span><<\/span>long<\/span> double<\/span>><\/span><\/span><\/span>(<\/span>vi)<\/span>;<\/span>
\n return<\/span> std::<\/span>sqrt<\/span>(<\/span>static_cast<\/span><<\/span>double<\/span>><\/span><\/span><\/span>(<\/span>s)<\/span>)<\/span>;<\/span>
\n }<\/span>;<\/span><\/p>\n

const<\/span> double<\/span> sc =<\/span> std::<\/span>max<\/span>(<\/span>atol,<\/span> rtol *<\/span> std::<\/span>max<\/span>(<\/span>norm2<\/span>(<\/span>y)<\/span>,<\/span> norm2<\/span>(<\/span>y_next)<\/span>)<\/span>)<\/span>;<\/span>
\n const<\/span> double<\/span> e2 =<\/span> norm2<\/span>(<\/span>err)<\/span>;<\/span>
\n return<\/span> (<\/span>sc ><\/span> 0.0<\/span>)<\/span> ?<\/span> (<\/span>e2 \/<\/span> sc)<\/span> :<\/span> e2;<\/span>
\n}<\/span><\/p>\n

\/\/ — C API wrappers (extern "C") —<\/span>
\nextern<\/span> "C"<\/span> {<\/span>
\ntypedef<\/span> void<\/span> (<\/span>*<\/span>ode45_deriv_cb)<\/span>(<\/span>double<\/span> t,<\/span> const<\/span> double<\/span>*<\/span> y,<\/span> double<\/span>*<\/span> dydt,<\/span> size_t n,<\/span> void<\/span>*<\/span> ctx)<\/span>;<\/span><\/p>\n

int<\/span> ode45_integrate_final<\/span>(<\/span>double<\/span> t0,<\/span>
\n const<\/span> double<\/span>*<\/span> y0,<\/span>
\n size_t n,<\/span>
\n double<\/span> t1,<\/span>
\n double<\/span> h_init,<\/span>
\n double<\/span> atol,<\/span>
\n double<\/span> rtol,<\/span>
\n int<\/span> norm_control,<\/span>
\n int<\/span> max_steps,<\/span>
\n ode45_deriv_cb cb,<\/span>
\n void<\/span>*<\/span> ctx,<\/span>
\n double<\/span>*<\/span> y_out)<\/span> {<\/span>
\n try<\/span> {<\/span>
\n State y0v<\/span>(<\/span>n)<\/span>;<\/span>
\n for<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> n;<\/span> ++<\/span>i)<\/span> y0v[<\/span>i]<\/span> =<\/span> y0[<\/span>i]<\/span>;<\/span>
\n DerivFunc f =<\/span> [<\/span>cb,<\/span> n,<\/span> ctx]<\/span>(<\/span>double<\/span> t,<\/span> const<\/span> State&<\/span> y,<\/span> State&<\/span> dydt)<\/span> {<\/span>
\n dydt.<\/span>resize<\/span>(<\/span>n)<\/span>;<\/span>
\n cb<\/span>(<\/span>t,<\/span> y.<\/span>data<\/span>(<\/span>)<\/span>,<\/span> dydt.<\/span>data<\/span>(<\/span>)<\/span>,<\/span> n,<\/span> ctx)<\/span>;<\/span>
\n }<\/span>;<\/span>
\n State y1 =<\/span> integrate_runge_kutta_dopri5<\/span>(<\/span>f,<\/span> t0,<\/span> y0v,<\/span> t1,<\/span> h_init,<\/span> atol,<\/span> rtol,<\/span> norm_control !=<\/span> 0<\/span>,<\/span> max_steps)<\/span>;<\/span>
\n for<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> n;<\/span> ++<\/span>i)<\/span> y_out[<\/span>i]<\/span> =<\/span> y1[<\/span>i]<\/span>;<\/span>
\n return<\/span> 0<\/span>;<\/span>
\n }<\/span> catch<\/span> (<\/span>.<\/span>.<\/span>.<\/span>)<\/span> {<\/span>
\n return<\/span> –<\/span>1<\/span>;<\/span>
\n }<\/span>
\n}<\/span><\/p>\n

int<\/span> ode45_integrate_times<\/span>(<\/span>double<\/span> t0,<\/span>
\n const<\/span> double<\/span>*<\/span> y0,<\/span>
\n size_t n,<\/span>
\n const<\/span> double<\/span>*<\/span> times,<\/span>
\n size_t m,<\/span>
\n double<\/span> h_init,<\/span>
\n double<\/span> atol,<\/span>
\n double<\/span> rtol,<\/span>
\n int<\/span> normControl,<\/span>
\n int<\/span> max_steps,<\/span>
\n ode45_deriv_cb cb,<\/span>
\n void<\/span>*<\/span> ctx,<\/span>
\n double<\/span>*<\/span> y_out)<\/span> {<\/span>
\n try<\/span> {<\/span>
\n State y0v<\/span>(<\/span>n)<\/span>;<\/span>
\n for<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> n;<\/span> ++<\/span>i)<\/span> y0v[<\/span>i]<\/span> =<\/span> y0[<\/span>i]<\/span>;<\/span>
\n std::<\/span>vector<<\/span>double<\/span>><\/span> tv<\/span>(<\/span>m)<\/span>;<\/span>
\n for<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> m;<\/span> ++<\/span>i)<\/span> tv[<\/span>i]<\/span> =<\/span> times[<\/span>i]<\/span>;<\/span>
\n DerivFunc f =<\/span> [<\/span>cb,<\/span> n,<\/span> ctx]<\/span>(<\/span>double<\/span> t,<\/span> const<\/span> State&<\/span> y,<\/span> State&<\/span> dydt)<\/span> {<\/span>
\n dydt.<\/span>resize<\/span>(<\/span>n)<\/span>;<\/span>
\n cb<\/span>(<\/span>t,<\/span> y.<\/span>data<\/span>(<\/span>)<\/span>,<\/span> dydt.<\/span>data<\/span>(<\/span>)<\/span>,<\/span> n,<\/span> ctx)<\/span>;<\/span>
\n }<\/span>;<\/span>
\n Trajectory traj =<\/span> integrate_runge_kutta_dopri5_at_times<\/span>(<\/span>f,<\/span> t0,<\/span> y0v,<\/span> tv,<\/span> h_init,<\/span> atol,<\/span> rtol,<\/span> normControl !=<\/span> 0<\/span>,<\/span> max_steps)<\/span>;<\/span>
\n if<\/span> (<\/span>traj.<\/span>t.<\/span>size<\/span>(<\/span>)<\/span> !=<\/span> m ||<\/span> traj.<\/span>y.<\/span>size<\/span>(<\/span>)<\/span> !=<\/span> m)<\/span> return<\/span> –<\/span>2<\/span>;<\/span>
\n \/\/ Row-major blocks: for each time i, copy n entries<\/span>
\n for<\/span> (<\/span>size_t i =<\/span> 0<\/span>;<\/span> i <<\/span> m;<\/span> ++<\/span>i)<\/span> {<\/span>
\n const<\/span> State&<\/span> yi =<\/span> traj.<\/span>y[<\/span>i]<\/span>;<\/span>
\n for<\/span> (<\/span>size_t j =<\/span> 0<\/span>;<\/span> j <<\/span> n;<\/span> ++<\/span>j)<\/span> {<\/span>
\n y_out[<\/span>i*<\/span>n +<\/span> j]<\/span> =<\/span> yi[<\/span>j]<\/span>;<\/span>
\n }<\/span>
\n }<\/span>
\n return<\/span> 0<\/span>;<\/span>
\n }<\/span> catch<\/span> (<\/span>.<\/span>.<\/span>.<\/span>)<\/span> {<\/span>
\n return<\/span> –<\/span>1<\/span>;<\/span>
\n }<\/span>
\n}<\/span>
\n}<\/span><\/p>\n

\/*
\ny(5) = 0.00673794718173955
\nexact = 0.006737946999085467
\nabs err = 1.98847773490845e-011
\nsamples (first 5):
\n0.00000000: 1.00000000
\n0.00927251: 0.99077034
\n0.01854503: 0.98162587
\n0.02781754: 0.97256580
\n0.03709006: 0.96358935
\nlength refine=4: 349, expected 349
\n*\/<\/span><\/p>\n","protected":false},"excerpt":{"rendered":"

#include #ifdef __cplusplus
\nextern \\”C\\” {
\n#endif\/\/ C callback for derivatives: fills dydt of length n given y of length n at time t.
\ntypedef void (*ode45_deriv_cb)(double t, const double* y, double* dydt, size_t n, void* ctx);\/\/ In<\/p>\n","protected":false},"author":2,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[55,1284,190,427],"topic":[],"class_list":["post-106719","post","type-post","status-publish","format-standard","hentry","category-server","tag-c","tag-matlab","tag-190","tag-427"],"yoast_head":"\nC++\u4ee3\u7801\u5b9e\u73b0MATLAB\u4e2d\u7684ode45\u51fd\u6570\u529f\u80fd - \u7f51\u7855\u4e92\u8054\u5e2e\u52a9\u4e2d\u5fc3<\/title>\n<meta name=\"robots\" content=\"index, follow, max-snippet:-1, max-image-preview:large, max-video-preview:-1\" \/>\n<link rel=\"canonical\" href=\"https:\/\/www.wsisp.com\/helps\/106719.html\" \/>\n<meta property=\"og:locale\" content=\"zh_CN\" \/>\n<meta property=\"og:type\" content=\"article\" \/>\n<meta property=\"og:title\" content=\"C++\u4ee3\u7801\u5b9e\u73b0MATLAB\u4e2d\u7684ode45\u51fd\u6570\u529f\u80fd - \u7f51\u7855\u4e92\u8054\u5e2e\u52a9\u4e2d\u5fc3\" \/>\n<meta property=\"og:description\" content=\"#include #ifdef __cplusplus extern "C" { #endif\/\/ C callback for derivatives: fills dydt of length n given y of length n at time t. typedef void (*ode45_deriv_cb)(double t, const double* y, double* dydt, size_t n, void* ctx);\/\/ In\" \/>\n<meta property=\"og:url\" content=\"https:\/\/www.wsisp.com\/helps\/106719.html\" \/>\n<meta property=\"og:site_name\" content=\"\u7f51\u7855\u4e92\u8054\u5e2e\u52a9\u4e2d\u5fc3\" \/>\n<meta property=\"article:published_time\" content=\"2026-09-17T22:40:32+00:00\" \/>\n<meta name=\"author\" content=\"admin\" \/>\n<meta name=\"twitter:card\" content=\"summary_large_image\" \/>\n<meta name=\"twitter:label1\" content=\"\u4f5c\u8005\" \/>\n\t<meta name=\"twitter:data1\" content=\"admin\" \/>\n\t<meta name=\"twitter:label2\" content=\"\u9884\u8ba1\u9605\u8bfb\u65f6\u95f4\" \/>\n\t<meta name=\"twitter:data2\" content=\"18 \u5206\" \/>\n<script type=\"application\/ld+json\" class=\"yoast-schema-graph\">{\"@context\":\"https:\/\/schema.org\",\"@graph\":[{\"@type\":\"WebPage\",\"@id\":\"https:\/\/www.wsisp.com\/helps\/106719.html\",\"url\":\"https:\/\/www.wsisp.com\/helps\/106719.html\",\"name\":\"C++\u4ee3\u7801\u5b9e\u73b0MATLAB\u4e2d\u7684ode45\u51fd\u6570\u529f\u80fd - \u7f51\u7855\u4e92\u8054\u5e2e\u52a9\u4e2d\u5fc3\",\"isPartOf\":{\"@id\":\"https:\/\/www.wsisp.com\/helps\/#website\"},\"datePublished\":\"2026-09-17T22:40:32+00:00\",\"dateModified\":\"2026-09-17T22:40:32+00:00\",\"author\":{\"@id\":\"https:\/\/www.wsisp.com\/helps\/#\/schema\/person\/358e386c577a3ab51c4493330a20ad41\"},\"breadcrumb\":{\"@id\":\"https:\/\/www.wsisp.com\/helps\/106719.html#breadcrumb\"},\"inLanguage\":\"zh-Hans\",\"potentialAction\":[{\"@type\":\"ReadAction\",\"target\":[\"https:\/\/www.wsisp.com\/helps\/106719.html\"]}]},{\"@type\":\"BreadcrumbList\",\"@id\":\"https:\/\/www.wsisp.com\/helps\/106719.html#breadcrumb\",\"itemListElement\":[{\"@type\":\"ListItem\",\"position\":1,\"name\":\"\u9996\u9875\",\"item\":\"https:\/\/www.wsisp.com\/helps\"},{\"@type\":\"ListItem\",\"position\":2,\"name\":\"C++\u4ee3\u7801\u5b9e\u73b0MATLAB\u4e2d\u7684ode45\u51fd\u6570\u529f\u80fd\"}]},{\"@type\":\"WebSite\",\"@id\":\"https:\/\/www.wsisp.com\/helps\/#website\",\"url\":\"https:\/\/www.wsisp.com\/helps\/\",\"name\":\"\u7f51\u7855\u4e92\u8054\u5e2e\u52a9\u4e2d\u5fc3\",\"description\":\"\u9999\u6e2f\u670d\u52a1\u5668_\u9999\u6e2f\u4e91\u670d\u52a1\u5668\u8d44\u8baf_\u670d\u52a1\u5668\u5e2e\u52a9\u6587\u6863_\u670d\u52a1\u5668\u6559\u7a0b\",\"potentialAction\":[{\"@type\":\"SearchAction\",\"target\":{\"@type\":\"EntryPoint\",\"urlTemplate\":\"https:\/\/www.wsisp.com\/helps\/?s={search_term_string}\"},\"query-input\":\"required name=search_term_string\"}],\"inLanguage\":\"zh-Hans\"},{\"@type\":\"Person\",\"@id\":\"https:\/\/www.wsisp.com\/helps\/#\/schema\/person\/358e386c577a3ab51c4493330a20ad41\",\"name\":\"admin\",\"image\":{\"@type\":\"ImageObject\",\"inLanguage\":\"zh-Hans\",\"@id\":\"https:\/\/www.wsisp.com\/helps\/#\/schema\/person\/image\/\",\"url\":\"https:\/\/gravatar.wp-china-yes.net\/avatar\/?s=96&d=mystery\",\"contentUrl\":\"https:\/\/gravatar.wp-china-yes.net\/avatar\/?s=96&d=mystery\",\"caption\":\"admin\"},\"sameAs\":[\"http:\/\/wp.wsisp.com\"],\"url\":\"https:\/\/www.wsisp.com\/helps\/author\/admin\"}]}<\/script>\n<!-- \/ Yoast SEO plugin. -->","yoast_head_json":{"title":"C++\u4ee3\u7801\u5b9e\u73b0MATLAB\u4e2d\u7684ode45\u51fd\u6570\u529f\u80fd - \u7f51\u7855\u4e92\u8054\u5e2e\u52a9\u4e2d\u5fc3","robots":{"index":"index","follow":"follow","max-snippet":"max-snippet:-1","max-image-preview":"max-image-preview:large","max-video-preview":"max-video-preview:-1"},"canonical":"https:\/\/www.wsisp.com\/helps\/106719.html","og_locale":"zh_CN","og_type":"article","og_title":"C++\u4ee3\u7801\u5b9e\u73b0MATLAB\u4e2d\u7684ode45\u51fd\u6570\u529f\u80fd - \u7f51\u7855\u4e92\u8054\u5e2e\u52a9\u4e2d\u5fc3","og_description":"#include #ifdef __cplusplus extern \"C\" { #endif\/\/ C callback for derivatives: fills dydt of length n given y of length n at time t. typedef void (*ode45_deriv_cb)(double t, const double* y, double* dydt, size_t n, void* ctx);\/\/ In","og_url":"https:\/\/www.wsisp.com\/helps\/106719.html","og_site_name":"\u7f51\u7855\u4e92\u8054\u5e2e\u52a9\u4e2d\u5fc3","article_published_time":"2026-09-17T22:40:32+00:00","author":"admin","twitter_card":"summary_large_image","twitter_misc":{"\u4f5c\u8005":"admin","\u9884\u8ba1\u9605\u8bfb\u65f6\u95f4":"18 \u5206"},"schema":{"@context":"https:\/\/schema.org","@graph":[{"@type":"WebPage","@id":"https:\/\/www.wsisp.com\/helps\/106719.html","url":"https:\/\/www.wsisp.com\/helps\/106719.html","name":"C++\u4ee3\u7801\u5b9e\u73b0MATLAB\u4e2d\u7684ode45\u51fd\u6570\u529f\u80fd - \u7f51\u7855\u4e92\u8054\u5e2e\u52a9\u4e2d\u5fc3","isPartOf":{"@id":"https:\/\/www.wsisp.com\/helps\/#website"},"datePublished":"2026-09-17T22:40:32+00:00","dateModified":"2026-09-17T22:40:32+00:00","author":{"@id":"https:\/\/www.wsisp.com\/helps\/#\/schema\/person\/358e386c577a3ab51c4493330a20ad41"},"breadcrumb":{"@id":"https:\/\/www.wsisp.com\/helps\/106719.html#breadcrumb"},"inLanguage":"zh-Hans","potentialAction":[{"@type":"ReadAction","target":["https:\/\/www.wsisp.com\/helps\/106719.html"]}]},{"@type":"BreadcrumbList","@id":"https:\/\/www.wsisp.com\/helps\/106719.html#breadcrumb","itemListElement":[{"@type":"ListItem","position":1,"name":"\u9996\u9875","item":"https:\/\/www.wsisp.com\/helps"},{"@type":"ListItem","position":2,"name":"C++\u4ee3\u7801\u5b9e\u73b0MATLAB\u4e2d\u7684ode45\u51fd\u6570\u529f\u80fd"}]},{"@type":"WebSite","@id":"https:\/\/www.wsisp.com\/helps\/#website","url":"https:\/\/www.wsisp.com\/helps\/","name":"\u7f51\u7855\u4e92\u8054\u5e2e\u52a9\u4e2d\u5fc3","description":"\u9999\u6e2f\u670d\u52a1\u5668_\u9999\u6e2f\u4e91\u670d\u52a1\u5668\u8d44\u8baf_\u670d\u52a1\u5668\u5e2e\u52a9\u6587\u6863_\u670d\u52a1\u5668\u6559\u7a0b","potentialAction":[{"@type":"SearchAction","target":{"@type":"EntryPoint","urlTemplate":"https:\/\/www.wsisp.com\/helps\/?s={search_term_string}"},"query-input":"required name=search_term_string"}],"inLanguage":"zh-Hans"},{"@type":"Person","@id":"https:\/\/www.wsisp.com\/helps\/#\/schema\/person\/358e386c577a3ab51c4493330a20ad41","name":"admin","image":{"@type":"ImageObject","inLanguage":"zh-Hans","@id":"https:\/\/www.wsisp.com\/helps\/#\/schema\/person\/image\/","url":"https:\/\/gravatar.wp-china-yes.net\/avatar\/?s=96&d=mystery","contentUrl":"https:\/\/gravatar.wp-china-yes.net\/avatar\/?s=96&d=mystery","caption":"admin"},"sameAs":["http:\/\/wp.wsisp.com"],"url":"https:\/\/www.wsisp.com\/helps\/author\/admin"}]}},"_links":{"self":[{"href":"https:\/\/www.wsisp.com\/helps\/wp-json\/wp\/v2\/posts\/106719","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.wsisp.com\/helps\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.wsisp.com\/helps\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.wsisp.com\/helps\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/www.wsisp.com\/helps\/wp-json\/wp\/v2\/comments?post=106719"}],"version-history":[{"count":0,"href":"https:\/\/www.wsisp.com\/helps\/wp-json\/wp\/v2\/posts\/106719\/revisions"}],"wp:attachment":[{"href":"https:\/\/www.wsisp.com\/helps\/wp-json\/wp\/v2\/media?parent=106719"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.wsisp.com\/helps\/wp-json\/wp\/v2\/categories?post=106719"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.wsisp.com\/helps\/wp-json\/wp\/v2\/tags?post=106719"},{"taxonomy":"topic","embeddable":true,"href":"https:\/\/www.wsisp.com\/helps\/wp-json\/wp\/v2\/topic?post=106719"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}