octave:2> realmax ans = 1.7977e+308 octave:3> realmax+realmax ans = Inf octave:4> (realmax+realmax)/2 ans = Inf octave:5> realmax+(realmax-realmax)/2 ans = 1.7977e+308 octave:6> odefun=@(t,y) [-2*y(2)*y(1);y(1)^2+y(3)^2-y(2)^2-1;-2*(y(1)+y(2))*y(3)]; octave:7> options.InitialStep=1e-2; octave:8> options.RelTol=1e-2; octave:9> options.AbsTol=1e-2; octave:10> y0=[0.5;2;10]; octave:11> tspan=[0,15]; octave:12> [tout,yout]=rk23(odefun,tspan,y0,options); warning: time step rejected warning: time step rejected warning: time step rejected warning: time step rejected warning: time step rejected warning: time step rejected warning: time step rejected warning: time step rejected octave:13> length(tout) ans = 167 octave:14> yout yout = 5.0000e-01 2.0000e+00 1.0000e+01 4.7607e-01 2.8776e+00 9.4288e+00 4.2295e-01 3.9913e+00 8.2494e+00 3.6625e-01 4.6341e+00 7.0507e+00 3.0704e-01 4.9512e+00 5.8386e+00 2.5256e-01 4.9925e+00 4.7504e+00 2.0709e-01 4.8550e+00 3.8592e+00 1.6639e-01 4.5918e+00 3.0743e+00 1.3033e-01 4.2307e+00 2.3883e+00 1.0082e-01 3.8235e+00 1.8340e+00 7.8176e-02 3.4203e+00 1.4132e+00 6.0370e-02 3.0265e+00 1.0853e+00 4.6343e-02 2.6478e+00 8.2899e-01 3.5327e-02 2.2876e+00 6.2911e-01 2.6736e-02 1.9473e+00 4.7418e-01 2.0103e-02 1.6269e+00 3.5520e-01 1.5049e-02 1.3248e+00 2.6500e-01 1.1381e-02 1.0477e+00 1.9978e-01 8.8058e-03 7.9405e-01 1.5414e-01 7.0245e-03 5.5209e-01 1.2263e-01 5.8844e-03 3.0814e-01 1.0246e-01 5.3819e-03 3.8886e-02 9.3432e-02 5.8918e-03 -3.0569e-01 1.0190e-01 8.7754e-03 -7.8822e-01 1.5099e-01 1.2350e-02 -1.1267e+00 2.1170e-01 2.0454e-02 -1.6481e+00 3.4862e-01 3.5854e-02 -2.3170e+00 6.0644e-01 5.2762e-02 -2.8499e+00 8.8664e-01 8.6237e-02 -3.6221e+00 1.4344e+00 1.3827e-01 -4.4377e+00 2.2708e+00 2.1082e-01 -5.1479e+00 3.4110e+00 3.0833e-01 -5.6117e+00 4.9001e+00 4.2735e-01 -5.6058e+00 6.6510e+00 5.4203e-01 -5.0325e+00 8.2578e+00 6.5278e-01 -3.7737e+00 9.7010e+00 7.2644e-01 -2.0827e+00 1.0530e+01 7.5527e-01 -1.2575e-01 1.0666e+01 7.3597e-01 1.6508e+00 1.0140e+01 6.8088e-01 3.0684e+00 9.1669e+00 6.0365e-01 4.0666e+00 7.9546e+00 5.2167e-01 4.6350e+00 6.7465e+00 4.3458e-01 4.8974e+00 5.5189e+00 3.5940e-01 4.8973e+00 4.4948e+00 2.9463e-01 4.7380e+00 3.6356e+00 2.3641e-01 4.4612e+00 2.8805e+00 1.8460e-01 4.0913e+00 2.2221e+00 1.4368e-01 3.6958e+00 1.7113e+00 1.1171e-01 3.3025e+00 1.3185e+00 8.6409e-02 2.9183e+00 1.0115e+00 6.6389e-02 2.5487e+00 7.7145e-01 5.0626e-02 2.1968e+00 5.8437e-01 3.8318e-02 1.8640e+00 4.3961e-01 2.8815e-02 1.5500e+00 3.2874e-01 2.1587e-02 1.2532e+00 2.4499e-01 1.6412e-02 9.8301e-01 1.8540e-01 1.2781e-02 7.3331e-01 1.4377e-01 1.0298e-02 4.9237e-01 1.1537e-01 8.7742e-03 2.4496e-01 9.7888e-02 8.2927e-03 -3.6471e-02 9.2082e-02 9.7353e-03 -4.1406e-01 1.0742e-01 1.2854e-02 -7.3610e-01 1.4107e-01 2.0866e-02 -1.2193e+00 2.2716e-01 3.6552e-02 -1.8162e+00 3.9389e-01 6.5228e-02 -2.5441e+00 6.9389e-01 9.5818e-02 -3.1067e+00 1.0084e+00 1.5561e-01 -3.9111e+00 1.6103e+00 2.4833e-01 -4.7541e+00 2.5161e+00 3.7906e-01 -5.4818e+00 3.7450e+00 5.5937e-01 -5.9341e+00 5.3569e+00 7.6913e-01 -5.8620e+00 7.1116e+00 9.7474e-01 -5.2046e+00 8.6873e+00 1.1778e+00 -3.8040e+00 1.0040e+01 1.3056e+00 -2.0744e+00 1.0658e+01 1.3560e+00 -2.7261e-02 1.0547e+01 1.3205e+00 1.6394e+00 9.8400e+00 1.2209e+00 2.9940e+00 8.7208e+00 1.0913e+00 3.8712e+00 7.5084e+00 9.4694e-01 4.3881e+00 6.2923e+00 7.8814e-01 4.6175e+00 5.0563e+00 6.5909e-01 4.5999e+00 4.1118e+00 5.4364e-01 4.4413e+00 3.3060e+00 4.3844e-01 4.1748e+00 2.6019e+00 3.4445e-01 3.8248e+00 1.9965e+00 2.7138e-01 3.4626e+00 1.5417e+00 2.1272e-01 3.0972e+00 1.1870e+00 1.6559e-01 2.7376e+00 9.0908e-01 1.2785e-01 2.3893e+00 6.9161e-01 9.7878e-02 2.0556e+00 5.2234e-01 7.4324e-02 1.7379e+00 3.9171e-01 5.6063e-02 1.4364e+00 2.9207e-01 4.2155e-02 1.1492e+00 2.1724e-01 3.2301e-02 8.8754e-01 1.6485e-01 2.5432e-02 6.4233e-01 1.2862e-01 2.0853e-02 4.0106e-01 1.0454e-01 1.8329e-02 1.4513e-01 9.1034e-02 1.8444e-02 -1.6111e-01 9.0596e-02 2.4439e-02 -5.9486e-01 1.1818e-01 3.3417e-02 -9.1897e-01 1.5970e-01 5.4739e-02 -1.4138e+00 2.5692e-01 9.5535e-02 -2.0390e+00 4.3807e-01 1.6964e-01 -2.8123e+00 7.5484e-01 3.0541e-01 -3.7747e+00 1.3050e+00 5.3107e-01 -4.8344e+00 2.1553e+00 8.7333e-01 -5.8495e+00 3.3275e+00 1.3889e+00 -6.6909e+00 4.8802e+00 1.8875e+00 -7.0324e+00 6.1706e+00 2.6043e+00 -6.9779e+00 7.6937e+00 3.4089e+00 -6.3210e+00 8.9321e+00 4.1646e+00 -5.1730e+00 9.5741e+00 4.8276e+00 -3.6342e+00 9.5489e+00 5.3150e+00 -1.8260e+00 8.7830e+00 5.4939e+00 -3.4847e-01 7.6726e+00 5.4453e+00 8.6982e-01 6.4162e+00 5.2195e+00 1.7963e+00 5.2055e+00 4.8155e+00 2.5389e+00 3.9931e+00 4.2915e+00 3.0281e+00 2.9509e+00 3.7285e+00 3.2817e+00 2.1456e+00 3.1961e+00 3.3577e+00 1.5668e+00 2.7152e+00 3.3177e+00 1.1524e+00 2.2683e+00 3.1946e+00 8.3904e-01 1.8600e+00 3.0068e+00 6.0289e-01 1.4942e+00 2.7686e+00 4.2647e-01 1.1740e+00 2.4926e+00 2.9645e-01 9.0084e-01 2.1901e+00 2.0220e-01 6.9047e-01 1.8964e+00 1.3969e-01 5.2813e-01 1.6142e+00 9.7323e-02 4.0330e-01 1.3440e+00 6.8250e-02 3.0833e-01 1.0846e+00 4.8217e-02 2.3746e-01 8.3307e-01 3.4463e-02 1.8649e-01 5.8350e-01 2.5172e-02 1.5392e-01 3.3370e-01 1.9341e-02 1.3884e-01 6.1753e-02 1.6153e-02 1.5043e-01 -2.8899e-01 1.5855e-02 2.2701e-01 -7.8594e-01 2.0796e-02 3.2534e-01 -1.1346e+00 2.6916e-02 5.5541e-01 -1.6686e+00 3.9072e-02 7.9321e-01 -2.0536e+00 4.9136e-02 1.2473e+00 -2.5776e+00 6.3620e-02 1.9217e+00 -3.0806e+00 7.7417e-02 2.8236e+00 -3.4461e+00 8.6545e-02 3.9996e+00 -3.5418e+00 8.7973e-02 5.1394e+00 -3.2482e+00 8.1140e-02 6.1943e+00 -2.5274e+00 6.8196e-02 6.9639e+00 -1.3936e+00 5.2042e-02 7.2459e+00 -1.0707e-01 3.7294e-02 7.0368e+00 1.2018e+00 2.4866e-02 6.4741e+00 2.2102e+00 1.6489e-02 5.6361e+00 2.9734e+00 1.0401e-02 4.8100e+00 3.3721e+00 6.8428e-03 3.9733e+00 3.5417e+00 4.4351e-03 3.2295e+00 3.5230e+00 2.9179e-03 2.6098e+00 3.3848e+00 1.9691e-03 2.0656e+00 3.1600e+00 1.3184e-03 1.5910e+00 2.8670e+00 8.6621e-04 1.2190e+00 2.5534e+00 5.7777e-04 9.3770e-01 2.2472e+00 3.9486e-04 7.1999e-01 1.9502e+00 2.7305e-04 5.5128e-01 1.6650e+00 1.9028e-04 4.2116e-01 1.3920e+00 1.3340e-04 3.2183e-01 1.1303e+00 9.4142e-05 2.4735e-01 8.7728e-01 6.7135e-05 1.9332e-01 6.2764e-01 4.8803e-05 1.5785e-01 3.7922e-01 3.7121e-05 1.3955e-01 1.1442e-01 3.0452e-05 1.4432e-01 -2.1241e-01 2.8732e-05 2.0470e-01 -6.8457e-01 3.5699e-05 3.3843e-01 -1.1779e+00 5.1213e-05 octave:15> E=(yout(:,1).^2+yout(:,2).^2+yout(:,3).^2+1)./(2*yout(:,1)) E = 105.2500 103.3563 100.6755 98.7313 97.2140 96.1249 95.3887 94.8462 94.4533 94.1930 94.0302 93.9288 93.8697 93.8413 93.8362 93.8497 93.8783 93.9104 93.9358 93.9485 93.9339 93.8587 93.6786 93.6805 93.7063 93.8262 93.9581 93.9227 93.8365 93.5458 92.9230 91.7811 89.9065 87.4595 84.0833 80.3521 76.3574 72.7497 69.6975 67.2381 65.4357 64.0063 63.0436 62.3708 61.8735 61.5119 61.2777 61.1270 61.0303 60.9708 60.9380 60.9251 60.9275 60.9419 60.9579 60.9703 60.9733 60.9552 60.8897 60.7626 60.7569 60.8372 60.9407 61.0036 60.9362 60.7755 60.3915 59.6464 58.3001 56.2508 53.6077 49.9447 46.1874 42.0620 38.7193 35.8358 33.6987 32.0741 30.7740 29.9656 29.3856 28.9565 28.6455 28.4472 28.3155 28.2280 28.1710 28.1356 28.1157 28.1072 28.1071 28.1092 28.1103 28.1063 28.0893 28.0447 27.9971 27.9967 28.0230 28.0444 28.0236 27.9052 27.5844 26.9378 25.7445 24.3955 22.2066 19.4135 16.4200 13.3290 10.3221 8.2067 6.6640 5.6104 4.8365 4.3451 4.0600 3.9023 3.8133 3.7594 3.7268 3.7075 3.6963 3.6904 3.6877 3.6867 3.6867 3.6874 3.6882 3.6889 3.6884 3.6855 3.6775 3.6775 3.6791 3.6857 3.6868 3.6896 3.6918 3.6931 3.6940 3.6941 3.6939 3.6934 3.6928 3.6921 3.6916 3.6911 3.6909 3.6910 3.6912 3.6915 3.6920 3.6928 3.6938 3.6947 3.6958 3.6969 3.6981 3.6994 3.7008 3.7019 3.7019 3.6997 3.6931 3.6896 3.6963 octave:16> plot(tout,E,'.') octave:17> stiff warning: axis: omitting non-positive data in log plot warning: axis: omitting non-positive data in log plot warning: axis: omitting non-positive data in log plot warning: axis: omitting non-positive data in log plot octave:18> stiff warning: axis: omitting non-positive data in log plot warning: axis: omitting non-positive data in log plot warning: axis: omitting non-positive data in log plot warning: axis: omitting non-positive data in log plot octave:19> close all octave:20> stiff warning: axis: omitting non-positive data in log plot warning: axis: omitting non-positive data in log plot warning: axis: omitting non-positive data in log plot warning: axis: omitting non-positive data in log plot octave:21> axis([1e2,1e4,1e-12,1]) octave:22> axis([1e2,1e4,1e-16,1]) octave:23> erree erree = Columns 1 through 5: 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 Columns 6 through 10: 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 Columns 11 through 15: 0.0000e+00 5.5720e+281 1.4933e+223 7.4185e+156 3.8443e+82 Columns 16 through 20: 1.0000e+00 1.3600e-18 1.3083e-18 1.2603e-18 1.2157e-18 Column 21: 1.1741e-18 octave:24> octave:24> octave:24> octave:24> octave:24> octave:24> octave:24> octave:24> octave:24> octave:24> octave:24> octave:24> octave:24> errei errei = Columns 1 through 6: 1.5166e-17 1.0984e-17 8.5285e-18 6.9347e-18 5.8254e-18 5.0127e-18 Columns 7 through 12: 4.3937e-18 3.9076e-18 3.5163e-18 3.1949e-18 2.9264e-18 2.6989e-18 Columns 13 through 18: 2.5038e-18 2.3347e-18 2.1867e-18 2.0562e-18 1.9403e-18 1.8366e-18 Columns 19 through 21: 1.7434e-18 1.6591e-18 1.5825e-18 octave:25> axis([1e2,1e4,1e-20,1e-4]) octave:26> axis([1e2,1e4,1e-20,1e-14]) octave:27> diary off