<div dir="ltr">Dear Scilabers<div><br></div><div>Just FYI, I had a scientific article published this weekend. You can see it  here:</div><div>















<span style="font-size:11pt;font-family:Calibri,sans-serif"><a href="http://www.aes.org/e-lib/browse.cfm?elib=19566" style="color:rgb(5,99,193);text-decoration:underline">http://www.aes.org/e-lib/browse.cfm?elib=19566</a></span>



<br></div><div><br></div><div>I'm not allowed to freely copy the article, but I hope nobody minds that I show the Scilab code included in the article. Here it is:</div><div><br></div><div><pre style="font-family:Monospaced;font-style:normal"><span style="color:rgb(0,0,0)">Qts</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(188,143,143)">0.5</span><span style="color:rgb(0,0,0)">;</span>
<span style="color:rgb(0,0,0)">alpha</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(188,143,143)">0</span><span style="color:rgb(0,0,0)">;</span> <span style="color:rgb(100,174,100);font-style:italic">// alpha = Vas/Vb</span>
           <span style="color:rgb(100,174,100);font-style:italic">// alpha = 0, open baffle</span>
<span style="color:rgb(0,0,0)">h</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(188,143,143)">0</span><span style="color:rgb(0,0,0)">;</span> <span style="color:rgb(100,174,100);font-style:italic">// bass reflex system tuning ratio</span>
       <span style="color:rgb(100,174,100);font-style:italic">// h = 0 for closed box systems</span>
<span style="color:rgb(176,24,19)">function</span> <span style="color:rgb(131,67,16);font-weight:bold">y</span><span style="color:rgb(92,92,92)">=</span><span style="color:rgb(0,0,0);text-decoration:underline">resp</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(131,67,16);font-weight:bold">s</span><span style="color:rgb(74,85,219)">)</span>
<span style="color:rgb(131,67,16);font-weight:bold">y</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(131,67,16);font-weight:bold">s</span><span style="color:rgb(92,92,92)">^</span><span style="color:rgb(188,143,143)">2</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(131,67,16);font-weight:bold">s</span><span style="color:rgb(92,92,92)">^</span><span style="color:rgb(188,143,143)">2</span><span style="color:rgb(92,92,92)">+</span><span style="color:rgb(131,67,16);font-weight:bold">s</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(0,0,0)">Qts</span><span style="color:rgb(92,92,92)">+</span><span style="color:rgb(188,143,143)">1</span><span style="color:rgb(92,92,92)">+</span><span style="color:rgb(0,0,0)">alpha</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
<span style="color:rgb(176,24,19)">endfunction</span> <span style="color:rgb(100,174,100);font-style:italic">// Eq (15)</span>

<span style="color:rgb(0,0,0)">N0</span>    <span style="color:rgb(92,92,92)">=</span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span> <span style="color:rgb(188,143,143)">12</span><span style="color:rgb(0,0,0)">;</span> <span style="color:rgb(100,174,100);font-style:italic">// resolution</span>
<span style="color:rgb(0,0,0)">steps</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(188,143,143)">8</span><span style="color:rgb(0,0,0)">;</span>  <span style="color:rgb(100,174,100);font-style:italic">// steps per unit time,</span>
<span style="color:rgb(0,0,0)">t_max</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(188,143,143)">5</span><span style="color:rgb(0,0,0)">;</span>  <span style="color:rgb(100,174,100);font-style:italic">// time (dimensionless)</span>
<span style="color:rgb(0,0,0)">t_step</span> <span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(188,143,143)">1</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(0,0,0)">steps</span><span style="color:rgb(0,0,0)">;</span>
<span style="color:rgb(0,0,0)">t</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(0,0,0)">t_step</span><span style="color:rgb(255,170,0)">:</span><span style="color:rgb(0,0,0)">t_step</span><span style="color:rgb(255,170,0)">:</span><span style="color:rgb(0,0,0)">t_max</span><span style="color:rgb(0,0,0)">;</span>

<span style="color:rgb(0,0,0)">mu_c</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(50,185,185)">max</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(74,85,219)">[</span><span style="color:rgb(188,143,143)">1</span> <span style="color:rgb(0,0,0)">h</span><span style="color:rgb(74,85,219)">]</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
<span style="color:rgb(0,0,0)">t_c</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(218,112,214)">%pi</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(0,0,0)">N0</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(188,143,143)">1</span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)"></span><span style="color:rgb(188,143,143)">2</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(92,92,92)"></span><span style="color:rgb(0,0,0)">mu_c</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
<span style="color:rgb(0,0,0)">nt</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(50,185,185)">length</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(0,0,0)">t</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>

<span style="color:rgb(50,185,185)">mprintf</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(188,143,143)">   t      -  Contour</span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
<span style="color:rgb(50,185,185)">mprintf</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(188,143,143)"> method  -  Exact\n</span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
<span style="color:rgb(160,32,240)">for</span> <span style="color:rgb(0,0,0)">i</span><span style="color:rgb(92,92,92)">=</span><span style="color:rgb(188,143,143)">1</span><span style="color:rgb(255,170,0)">:</span><span style="color:rgb(0,0,0)">nt</span> <span style="color:rgb(160,32,240)">do</span>
    <span style="color:rgb(0,0,0)">ti</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(0,0,0)">t</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)">i</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
    <span style="color:rgb(100,174,100);font-style:italic">// Selection defined in</span>
    <span style="color:rgb(100,174,100);font-style:italic">// Eqs (31) and (32)</span>
    <span style="color:rgb(160,32,240)">if</span> <span style="color:rgb(0,0,0)">ti</span> <span style="color:rgb(92,92,92)"><</span> <span style="color:rgb(0,0,0)">t_c</span> <span style="color:rgb(160,32,240)">then</span>
        <span style="color:rgb(0,0,0)">mu</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(218,112,214)">%pi</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(0,0,0)">N</span><span style="color:rgb(0,0,0)">0</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(188,143,143)">1</span><span style="color:rgb(188,143,143)">2</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(0,0,0)">ti</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
        <span style="color:rgb(0,0,0)">N</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(0,0,0)">N0</span><span style="color:rgb(0,0,0)">;</span>
        <span style="color:rgb(0,0,0)">d</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(188,143,143)">3</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(0,0,0)">N</span><span style="color:rgb(0,0,0)">;</span>
    <span style="color:rgb(160,32,240)">else</span>
        <span style="color:rgb(0,0,0)">mu</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(0,0,0)">mu_c</span><span style="color:rgb(0,0,0)">;</span>
        <span style="color:rgb(0,0,0)">N</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(50,185,185)">ceil</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(0,0,0)">N0</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(0,0,0)">ti</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(0,0,0)">t_c</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
        <span style="color:rgb(0,0,0)">d</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(188,143,143)">3</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(0,0,0)">N</span><span style="color:rgb(0,0,0)">;</span>
    <span style="color:rgb(160,32,240)">end</span>

    <span style="color:rgb(100,174,100);font-style:italic">// Perform sum in Eq (24)</span>
    <span style="color:rgb(0,0,0)">xsk</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(188,143,143)">0.0</span><span style="color:rgb(0,0,0)">;</span>
    <span style="color:rgb(160,32,240)">fo</span><span style="color:rgb(160,32,240)">r</span> <span style="color:rgb(0,0,0)">k</span><span style="color:rgb(92,92,92)">=</span><span style="color:rgb(92,92,92)">-</span><span style="color:rgb(0,0,0)">N</span><span style="color:rgb(255,170,0)">:</span><span style="color:rgb(0,0,0)"></span><span style="color:rgb(0,0,0)">N</span> <span style="color:rgb(160,32,240)">do</span>
        <span style="color:rgb(0,0,0)">u</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(0,0,0)">k</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(0,0,0)">d</span><span style="color:rgb(0,0,0)">;</span>
        <span style="color:rgb(0,0,0)">s</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(0,0,0)">mu</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(218,112,214)">%i</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(0,0,0)">u</span><span style="color:rgb(92,92,92)">+</span><span style="color:rgb(188,143,143)">1</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">^</span><span style="color:rgb(188,143,143)">2</span><span style="color:rgb(0,0,0)">;</span>
        <span style="color:rgb(0,0,0)">dsdu</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(188,143,143)">2</span><span style="color:rgb(188,143,143)">.0</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(0,0,0)">mu</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(218,112,214)">%i</span><span style="color:rgb(92,92,92)">-</span><span style="color:rgb(0,0,0)">u</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
        <span style="color:rgb(0,0,0)">xsk</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(0,0,0)">xs</span><span style="color:rgb(0,0,0)">k</span><span style="color:rgb(92,92,92)">+</span><span style="color:rgb(50,185,185)">imag</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(50,185,185)">exp</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(0,0,0)">s</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(0,0,0)">ti</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">*</span> <span style="color:rgb(255,170,0)">..</span>
            <span style="color:rgb(0,0,0);text-decoration:underline">resp</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(0,0,0)">s</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(0,0,0)">s</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(0,0,0)">dsdu</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
    <span style="color:rgb(160,32,240)">end</span>
    <span style="color:rgb(0,0,0)">xsk</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(0,0,0)">xsk</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(0,0,0)">d</span><span style="color:rgb(92,92,92)">/</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(188,143,143)">2</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(218,112,214)">%pi</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
    <span style="color:rgb(0,0,0)">xs</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(0,0,0)">i</span><span style="color:rgb(74,85,219)">)</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(0,0,0)">xsk</span><span style="color:rgb(0,0,0)">;</span> <span style="color:rgb(100,174,100);font-style:italic">// Collect x-values</span>
    <span style="color:rgb(0,0,0)">exact</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(0,0,0)">i</span><span style="color:rgb(74,85,219)">)</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(50,185,185)">exp</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(92,92,92)">-</span><span style="color:rgb(0,0,0)">ti</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(92,92,92)">*</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(188,143,143)">1</span><span style="color:rgb(92,92,92)">-</span><span style="color:rgb(0,0,0)">ti</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
    <span style="color:rgb(100,174,100);font-style:italic">// Collect exact result, compare</span>
    <span style="color:rgb(50,185,185)">mprintf</span><span style="color:rgb(74,85,219)"></span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(188,143,143)">%</span><span style="color:rgb(188,143,143)">5.4e </span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(0,0,0)">ti</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
    <span style="color:rgb(50,185,185)">mprintf</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(188,143,143)">%+12.11e </span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(0,0,0)">xsk</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
    <span style="color:rgb(50,185,185)">mprintf</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(188,143,143)">%+12.11e\n</span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(0,0,0)">exact</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(0,0,0)">i</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
<span style="color:rgb(160,32,240)">end</span>

<span style="color:rgb(0,0,0)">g_tr</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(0,0,0)">scf</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span> <span style="color:rgb(100,174,100);font-style:italic">// Graph Time Response</span>
<span style="color:rgb(0,0,0)">plot</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(74,85,219)">[</span><span style="color:rgb(188,143,143)">0</span> <span style="color:rgb(0,0,0)">t</span><span style="color:rgb(74,85,219)">]</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(74,85,219)">[</span><span style="color:rgb(188,143,143)">1</span><span style="color:rgb(0,0,0)">;</span> <span style="color:rgb(0,0,0)">xs</span><span style="color:rgb(74,85,219)">]</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(188,143,143)">-xr</span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
<span style="color:rgb(50,185,185)">xgrid</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(50,185,185)">color</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(188,143,143)">grey70</span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
<span style="color:rgb(100,174,100);font-style:italic">// Plot contour method numerical error</span>
<span style="color:rgb(0,0,0)">err</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(50,185,185)">abs</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(0,0,0)">xs</span> <span style="color:rgb(92,92,92)">-</span> <span style="color:rgb(0,0,0)">exact</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
<span style="color:rgb(0,0,0)">g_e</span> <span style="color:rgb(92,92,92)">=</span> <span style="color:rgb(0,0,0)">scf</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span> <span style="color:rgb(100,174,100);font-style:italic">// Graph Error</span>
<span style="color:rgb(0,0,0)">plot</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(74,85,219)">[</span><span style="color:rgb(188,143,143)">0</span> <span style="color:rgb(0,0,0)">t</span><span style="color:rgb(74,85,219)">]</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(74,85,219)">[</span><span style="color:rgb(188,143,143)">0</span><span style="color:rgb(0,0,0)">;</span> <span style="color:rgb(0,0,0)">err</span><span style="color:rgb(74,85,219)">]</span><span style="color:rgb(0,0,0)">,</span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(188,143,143)">-r</span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span>
<span style="color:rgb(50,185,185)">xgrid</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(50,185,185)">color</span><span style="color:rgb(74,85,219)">(</span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(188,143,143)">grey70</span><span style="color:rgb(188,143,143)">"</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(74,85,219)">)</span><span style="color:rgb(0,0,0)">;</span></pre><br></div><div><br></div><div>Best regards,</div><div>Claus</div></div>