1function varargout = ode2345(svname,ode,tspan,y0,options,varargin)
  2%MATLAB Code Generation Private Function
  3
  4%   Copyright 1984-2015 The MathWorks, Inc.
  5%#codegen
  6
  7%--------------------------------------------------------------------------
  8% Validate output mode
  9GraphicalOutput = nargout == 0; % odeXX(...)
 10StructOutput = nargout == 1; % sol = odeXX(...)
 11ArrayOutput  = ~GraphicalOutput && ~StructOutput; % [t,y, ...] = odeXX(...)
 12% For code generation,{0} requires at least {1,number,integer} outputs.
 13% Options corresponding to fewer outputs are not supported.
 14eml_invariant(ArrayOutput, ...
 15    'Coder:toolbox:RequiresNOutputs',svname,2);
 16% GrowMainArrays = StructOutput || numel(tspan) <= 2;
 17% Variable sizing support is required even when tspan has more than two
 18% elements because terminating because of minimum step size can result in a
 19% shorter array.
 20eml_invariant(eml_option('VariableSizing'), ...
 21    'Coder:toolbox:RequiresVariableSizing',svname);
 22
 23%--------------------------------------------------------------------------
 24% Validate inputs
 25eml_invariant(nargin >= 5,'MATLAB:minrhs');
 26coder.internal.prefer_const(svname,ode,tspan,options);
 27eml_invariant(eml_is_const(svname) && ...
 28    (strcmp(svname,'ode23') || strcmp(svname,'ode45')), ...
 29    'Coder:toolbox:UnsupportedSolver');
 30ODE23 = coder.const(strcmp(svname,'ode23'));
 31eml_invariant((ischar(ode) && eml_is_const(ode)) || ...
 32    isa(ode,'function_handle'), ...
 33    'Coder:toolbox:UnsupportedODEFunction');
 34if ischar(ode)
 35    ODEFunction = str2func(ode);
 36else
 37    ODEFunction = ode;
 38end
 39eml_invariant(isreal(tspan), ...
 40    'Coder:toolbox:TspanMustBeReal');
 41coder.internal.errorIf(isempty(tspan) || isempty(y0), ...
 42    'MATLAB:odearguments:TspanOrY0NotSupplied',svname);
 43coder.internal.errorIf(length(tspan) < 2, ...
 44    'MATLAB:odearguments:SizeTspan',svname);
 45t0 = tspan(1);
 46tfinal = tspan(end);
 47coder.internal.errorIf(t0 == tfinal, ...
 48    'MATLAB:odearguments:TspanEndpointsNotDistinct');
 49eml_invariant(ismonotonic(tspan), ...
 50    'MATLAB:odearguments:TspanNotMonotonic');
 51eml_invariant(isvector(y0), ...
 52    'Coder:toolbox:y0MustBeAVector');
 53
 54%--------------------------------------------------------------------------
 55% Validate mass matrix option
 56Moption = odeget(options,'Mass','','fast');
 57if ~isempty(Moption) % non-trivial mass matrix
 58    Msingular = odeget(options,'MassSingular','no','fast');
 59    coder.internal.errorIf(strcmp(Msingular,'yes'), ...
 60        make_error_id(svname,'MassSingularYes'));
 61    if strcmp(Msingular,'maybe')
 62        coder.internal.warning(make_error_id(svname,'MassSingularAssumedNo'));
 63    end
 64end
 65
 66%--------------------------------------------------------------------------
 67% Miscellaneous initializations
 68ZERO = coder.internal.indexInt(0);
 69ONE = coder.internal.indexInt(1);
 70nsteps  = ZERO;
 71nfailed = ZERO;
 72nfevals = ZERO;
 73neq = coder.internal.indexInt(numel(y0));
 74
 75%--------------------------------------------------------------------------
 76% Make first function evaluation and determine output type.
 77% Non-negative solution components
 78idxNonNegative = getIdxNonNegative(options,y0);
 79AnyNonNegative = ~isempty(idxNonNegative);
 80% Before the first call to ODEFunction, make sure that at least the first
 81% two inputs to ODEFunction are the same class.
 82eml_invariant(isfloat(y0) && isa(y0,class(t0)), ...
 83    'Coder:toolbox:InconsistentTypes');
 84f0 = callODEFunctionNSM(ODEFunction,t0,y0(:),options,varargin{:});
 85nfevals = nfevals + 1;
 86datatype = coder.internal.scalarEg(tspan,y0,f0);
 87realtype = real(datatype);
 88realnull = zeros(0,0,'like',realtype);
 89datanull = zeros(0,0,'like',datatype);
 90outcls = class(datatype);
 91% The tspan and y0 inputs and the output of the ode function must have the
 92% same type, all single or all double.
 93eml_invariant(isfloat(datatype) && ...
 94    isa(tspan,outcls) && isa(y0,outcls) && isa(f0,outcls), ...
 95    'Coder:toolbox:InconsistentTypes');
 96[m,n] = size(f0);
 97% The output of the ODE function must be a column vector.
 98eml_invariant(eml_is_const(n) && n == 1, ...
 99    'Coder:toolbox:FOMustReturnCol');
100% The size of y0 and the size of the output of the ODE function must match.
101eml_invariant(m == neq, ...
102    'Coder:toolbox:SizeIC')
103
104%--------------------------------------------------------------------------
105% Set tolerances and threshold
106rtol = cast(odeget(options,'RelTol',1e-3,'fast'),'like',realtype);
107eml_invariant(isscalar(rtol) && rtol > 0, ...
108    'MATLAB:odearguments:RelTolNotPosScalar');
109minrtol = 100*eps(outcls);
110if rtol < minrtol
111    rtol = minrtol;
112    coder.internal.warning('MATLAB:odearguments:RelTolIncrease', ...
113        coder.internal.num2str(rtol));
114end
115atol = cast(odeget(options,'AbsTol',1e-6,'fast'),'like',realtype);
116eml_invariant(ispositive(atol), ...
117    'MATLAB:odearguments:AbsTolNotPos');
118normcontrol = strcmp(odeget(options,'NormControl','off','fast'),'on');
119if normcontrol
120    eml_invariant(isscalar(atol), ...
121        'MATLAB:odearguments:NonScalarAbsTol');
122    normy = norm(y0(:));
123else
124    % AbsTol must be a scalar or a vector of length {0,number,integer}.
125    eml_invariant(isscalar(atol) || numel(atol) == neq, ...
126        'Coder:toolbox:SizeAbsTol',neq);
127    normy = zeros('like',realtype);
128end
129threshold = atol(:)/rtol;
130
131%--------------------------------------------------------------------------
132% Initialize method parameters.
133if ODE23
134    nstages = coder.internal.indexInt(4);
135    pow = cast(1/3,'like',realtype);
136    A = cast([1/2, 3/4, 1],'like',realtype);
137    % B = [
138    %     1/2 0   2/9
139    %     0   3/4 1/3
140    %     0   0   4/9
141    %     0   0   0
142    %     ];
143    % Coefficients of the explicit formula in packed form.
144    B = cast( ...
145        [1/2; ...
146        0; 3/4; ...
147        2/9; 1/3; 4/9], ...
148        outcls);
149    E = cast([-5/72; 1/12; 1/9; -1/8],'like',realtype);
150    ntrpfun = @ntrp23;
151    MinStepFactor = cast(0.5,'like',realtype);
152    DefaultRefine = 1;
153else % ODE45
154    nstages = coder.internal.indexInt(7);
155    pow = cast(1/5,'like',realtype);
156    A = cast([1/5,3/10,4/5,8/9,1,1],'like',realtype);
157    % B = [
158    %     1/5  3/40  44/45  19372/6561  9017/3168     35/384
159    %     0    9/40 -56/15 -25360/2187  -355/33        0
160    %     0    0     32/9   64448/6561 46732/5247    500/1113
161    %     0    0      0      -212/729     49/176     125/192
162    %     0    0      0         0      -5103/18656 -2187/6784
163    %     0    0      0         0          0          11/84
164    %     0    0      0         0          0           0
165    %     ];
166    % Coefficients of the explicit formula in packed form.
167    B = cast( ...
168        [1/5; ...
169        3/40; 9/40; ...
170        44/45; -56/15; 32/9; ...
171        19372/6561; -25360/2187; 64448/6561; -212/729; ...
172        9017/3168; -355/33; 46732/5247; 49/176; -5103/18656; ...
173        35/384; 0; 500/1113; 125/192; -2187/6784; 11/84], ...
174        'like',realtype);
175    E = cast( ...
176        [71/57600; 0; -71/16695; 71/1920; -17253/339200; 22/525; -1/40],...
177        'like',realtype);
178    ntrpfun = @ntrp45;
179    MinStepFactor = cast(0.1,'like',realtype);
180    DefaultRefine = 4;
181end
182
183%--------------------------------------------------------------------------
184% Process event function option
185eventOption = odeget(options,'Events','','fast');
186HasEvents = ~isempty(eventOption);
187if HasEvents
188    % You must enable dynamic memory allocation to use an event function.
189    eml_invariant(~strcmpi(eml_option('UseMalloc'),'Off'), ...
190        'Coder:toolbox:EventsRequireDynamicMemoryAllocation',svname);
191    if ischar(eventOption)
192        EventFunction = str2func(eventOption);
193    else
194        EventFunction = eventOption;
195    end
196    valt = EventFunction(t0,y0(:),varargin{:});
197end
198
199%--------------------------------------------------------------------------
200% Prepare output arrays
201% Allocate memory if we're generating output.
202teout = zeros(1,0,'like',realtype);
203yeout = zeros(neq,0,'like',datatype);
204ieout = zeros(1,0,coder.internal.indexIntClass);
205ntspan = numel(tspan);
206refine = coder.internal.indexInt( ...
207    max(1,odeget(options,'Refine',DefaultRefine,'fast')));
208if eml_is_const(ntspan) && ntspan == 2
209    coder.varsize('tout',[],[0,1]);
210    coder.varsize('yout');
211end
212if eml_is_const(size(tspan)) && ntspan == 2 % alloc in chunks
213    % Note that chunk is guaranteed to be greater than refine.
214    chunk = coder.internal.indexInt( ...
215        min(max(100,50*refine),refine+floor((2^13)/neq)));
216    tout = zeros(1,chunk,'like',realtype);
217    yout = zeros(neq,chunk,'like',datatype);
218    if refine <= 1
219        outputAt = SOLVERSTEPS; % computed points, no refinement
220    else
221        outputAt = REFINEDSTEPS; % computed points, with refinement
222        S = cast(1:refine-1,'like',realtype)/cast(refine,'like',realtype);
223    end
224else % output only at tspan points
225    outputAt = REQUESTEDPOINTS; % output only at tspan points
226    if ntspan == 2
227        % The tspan input is variable-size and became a 2-element array at
228        % run-time. Results will not match MATLAB in this case.
229        coder.internal.warning('Coder:toolbox:VarsizeTspanHasLengthTwo');
230    end
231    tout = zeros(1,ntspan,'like',realtype);
232    yout = zeros(neq,ntspan,'like',datatype);
233end
234nout = ONE;
235tout(nout) = t0;
236yout(:,nout) = y0(:);
237
238%--------------------------------------------------------------------------
239% Initialize output function
240outputFunctionOption = odeget(options,'OutputFcn','','fast');
241HasOutputFunction = ~isempty(outputFunctionOption);
242if HasOutputFunction
243    if ischar(outputFunctionOption)
244        OutputFunction = str2func(outputFunctionOption);
245    else
246        OutputFunction = outputFunctionOption;
247    end
248    outputidx = odeget(options,'OutputSel',ONE:neq,'fast');
249    eml_invariant(isValidIdxVector(outputidx,neq), ...
250        'Coder:toolbox:InvalidOutputSel');
251    noutputidx = coder.internal.indexInt(numel(outputidx));
252    OutputFunction([t0,tfinal],y0(outputidx(:)),'init',varargin{:});
253end
254
255%--------------------------------------------------------------------------
256% Printstats not supported
257PrintStats = strcmp(odeget(options,'Stats','off','fast'),'on');
258% Printing of output statistics is not supported for code generation.
259eml_invariant(~PrintStats, ...
260    'Coder:toolbox:NoPrintStats');
261
262%--------------------------------------------------------------------------
263% Compute/validates step size limits.
264% By default, hmax is 1/10 of the interval.
265hmax = min( ...
266    abs(tfinal - t0), ...
267    abs(odeget(options,'MaxStep',0.1*(tfinal - t0),'fast')));
268eml_invariant(hmax > 0, ...
269    'MATLAB:odearguments:MaxStepLEzero');
270htry = abs(odeget(options,'InitialStep',[],'fast'));
271eml_invariant(isempty(htry) || htry > 0, ...
272    'MATLAB:odearguments:InitialStepLEzero');
273hmin = 16*eps(t0);
274if isempty(htry)
275    % Compute an initial step size h using y'(t).
276    htspan = abs(tspan(2) - tspan(1));
277    absh = min(hmax,htspan);
278    if normcontrol
279        rh = (norm(f0)/max(normy,threshold))/(0.8*rtol^pow);
280    else
281        rh = norm(f0./max(abs(y0(:)),threshold),'inf')/(0.8*rtol^pow);
282    end
283    if absh*rh > 1
284        absh = 1/rh;
285    end
286    absh = max(absh,hmin);
287else
288    absh = min(hmax,max(hmin,htry));
289end
290
291%--------------------------------------------------------------------------
292% Final variable initializations
293t = t0;
294y = cast(y0(:),'like',datatype); % y0 might be real while y can be complex.
295f = zeros(neq,nstages,'like',y);
296f(:,1) = f0;
297tdir = sign(tfinal - t0);
298next = coder.internal.indexInt(2); % next entry in tspan
299TerminalEvent = false; % Compiler currently requires this.
300
301%--------------------------------------------------------------------------
302% THE MAIN LOOP
303MinStepExit = false;
304Done = false;
305while ~Done
306    % By default, hmin is a small number such that t+hmin is only slightly
307    % different than t.  It might be 0 if t is 0.
308    hmin = 16*eps(t);
309    absh = min(hmax,max(hmin,absh)); % couldn't limit absh until new hmin
310    h = tdir*absh;
311    % Stretch the step if within 10% of tfinal-t.
312    if 1.1*absh >= abs(tfinal - t)
313        h = tfinal - t;
314        absh = abs(h);
315        Done = true;
316    end
317    % LOOP FOR ADVANCING ONE STEP.
318    NoFailedAttempts = true;
319    NNresetLastStage = false;
320    while true
321        Bcolidx = ONE;
322        for j = 1:nstages-2
323            % ystage = y + f*(h*Bcol);
324            Bcolidx = Bcolidx + j - 1;
325            ystage = coder.internal.blas.xgemv( ...
326                'N',neq,j,h,f,1,neq,B,Bcolidx,1,1,y,1,1);
327            f(:,j+1) = callODEFunctionNSM( ...
328                ODEFunction,t+h*A(j),ystage,options,varargin{:});
329        end
330        tnew = t + h*A(nstages-1);
331        if Done
332            tnew = tfinal; % Hit end point exactly.
333        end
334        % ynew = y + f*(h*B6);
335        Bcolidx = Bcolidx + nstages - 2;
336        ynew = coder.internal.blas.xgemv( ...
337            'N',neq,nstages-1,h,f,1,neq,B,Bcolidx,1,1,y,1,1);
338        f(:,nstages) = callODEFunctionNSM( ...
339            ODEFunction,tnew,ynew,options,varargin{:});
340        nfevals = nfevals + (nstages - 1);
341        h = tnew - t; % Purify h.
342        % Estimate the error.
343        if normcontrol
344            normynew = norm(ynew);
345            errwt = max3(normy,normynew,threshold);
346            err = absh*(norm(f*E)/errwt);
347        else
348            err = absh*norm((f*E)./maxAbsThresh(y,ynew,threshold),'inf');
349            errwt = threshold;
350        end
351        NNrejectStep = false;
352        if AnyNonNegative && err <= rtol
353            errNN = normNegative(ynew,errwt,idxNonNegative,normcontrol);
354            if errNN > rtol
355                err = errNN;
356                NNrejectStep = true;
357            end
358        end
359        % Accept the solution only if the weighted error is no more than
360        % the tolerance rtol.  Estimate an h that will yield an error of
361        % rtol on the next step or the next try at taking this step, as the
362        % case may be, and use 0.8 of this value to avoid failures.
363        if err > rtol % Failed step
364            nfailed = nfailed + 1;
365            if absh <= hmin
366                coder.internal.warning( ...
367                    make_error_id(svname,'IntegrationTolNotMet'), ...
368                    coder.internal.num2str(t), ...
369                    coder.internal.num2str(hmin));
370                MinStepExit = true;
371                break
372            else
373                if NoFailedAttempts
374                    NoFailedAttempts = false;
375                    if NNrejectStep
376                        absh = max(hmin,0.5*absh);
377                    else
378                        absh = max(hmin, ...
379                            absh*max(MinStepFactor,0.8*(rtol/err)^pow));
380                    end
381                else
382                    absh = max(hmin,0.5*absh);
383                end
384                h = tdir*absh;
385                Done = false;
386            end
387        else % Successful step
388            NNresetLastStage = false;
389            if AnyNonNegative
390                for k = 1:numel(idxNonNegative)
391                    if ynew(idxNonNegative(k)) < 0
392                        ynew(idxNonNegative(k)) = 0;
393                        NNresetLastStage = true;
394                    end
395                end
396                if normcontrol
397                    normynew = norm(ynew);
398                end
399            end
400            break
401        end
402    end
403    if MinStepExit
404        break
405    end
406    nsteps = nsteps + 1;
407    if HasEvents
408        [ne,te,ye,ie,valt,TerminalEvent(1)] = ...
409            odezero(ntrpfun,EventFunction,valt,t,y,tnew,ynew,t0,h,f, ...
410            idxNonNegative,varargin{:});
411        if ne > 0
412            if StructOutput || nargout > 2
413                teout = [teout,te]; %#ok<AGROW>
414                yeout = [yeout,ye]; %#ok<AGROW>
415                ieout = [ieout,double(ie)]; %#ok<AGROW>
416            end
417            if TerminalEvent % Stop on a terminal event.
418                % Adjust the interpolation data to [t te(end)]. Update the
419                % derivatives using the interpolating polynomial.
420                taux = t + (te(end) - t)*A;
421                [~,f(:,2:nstages)] = ntrpfun(taux,t,y, ...
422                    realnull,datanull,h,f,idxNonNegative);
423                tnew = te(end);
424                ynew = ye(:,end);
425                h = tnew - t;
426                Done = true;
427            end
428        end
429    end
430    switch outputAt
431        case SOLVERSTEPS % computed points, no refinement
432            nout = nout + 1;
433            if nout > length(tout)
434                tout = [tout,zeros(1,chunk,'like',t)];  %#ok<AGROW>
435                yout = [yout,zeros(neq,chunk,'like',y)]; %#ok<AGROW>
436            end
437            tout(nout) = tnew;
438            for j = 1:neq
439                yout(j,nout) = ynew(j);
440            end
441            noutnew = ONE;
442            outidx = nout;
443        case REFINEDSTEPS % computed points, with refinement
444            noutnew = refine;
445            outidx = nout + 1;
446            tref = t + (tnew - t)*S;
447            toutnew = [tref,tnew];
448            youtnew = [ntrpfun(tref,t,y, ...
449                realnull,datanull,h,f,idxNonNegative),ynew];
450            nout = nout + noutnew;
451            if nout > length(tout)
452                % requires chunk >= refine
453                tout = [tout,zeros(1,chunk,'like',t)];  %#ok<AGROW>
454                yout = [yout,zeros(neq,chunk,'like',y)]; %#ok<AGROW>
455            end
456            % tout(firstout:nout) = tout_new;
457            % yout(:,firstout:nout) = yout_new;
458            for k = 1:noutnew
459                tout(k - 1 + outidx) = toutnew(k);
460                for j = 1:neq
461                    yout(j,k - 1 + outidx) = youtnew(j,k);
462                end
463            end
464        case REQUESTEDPOINTS % output only at tspan points
465            nnxt = NextNext(next,ntspan,tspan,tdir,tnew);
466            noutnew = nnxt - next;
467            outidx = next;
468            tnewused = false;
469            if noutnew > 0
470                i2 = nnxt - 1;
471                for k = outidx:i2-1
472                    tout(k) = tspan(k);
473                    yout(:,k) = ntrpfun(tspan(k),t,y, ...
474                        realnull,datanull,h,f,idxNonNegative);
475                end
476                tout(i2) = tspan(i2);
477                if tspan(i2) == tnew
478                    yout(:,i2) = ynew;
479                    tnewused = true;
480                else
481                    yout(:,i2) = ntrpfun(tspan(i2),t,y, ...
482                        realnull,datanull,h,f,idxNonNegative);
483                end
484                nout = nout + noutnew;
485                next = nnxt;
486            end
487            if HasEvents && TerminalEvent && ~tnewused && next <= ntspan
488                % Output tstop, ystop
489                tout(next) = tnew;
490                yout(:,next) = ynew;
491                nout = nout + 1;
492            end
493    end
494    if HasOutputFunction
495        toutnew = zeros(1,noutnew,'like',realtype);
496        youtnew = zeros(noutputidx,noutnew,'like',datatype);
497        for k = 1:noutnew
498            toutnew(k) = tout(k - 1 + outidx);
499            for j = 1:noutputidx
500                youtnew(j,k) = yout(outputidx(j),k - 1 + outidx);
501            end
502        end
503        TerminalEvent(1) = OutputFunction( ...
504            toutnew,youtnew,'',varargin{:});
505        if TerminalEvent
506            Done = true;
507        end
508    end
509    if Done
510        break
511    end
512    % If there were no failures compute a new h.
513    if NoFailedAttempts
514        % Note that absh may shrink by 0.8,and that err may be 0.
515        temp = 1.25*(err/rtol)^pow;
516        if temp > 0.2
517            absh = absh/temp;
518        else
519            absh = 5*absh;
520        end
521    end
522    % Advance the integration one step.
523    t = tnew;
524    y = ynew;
525    if normcontrol
526        normy = normynew;
527    end
528    if NNresetLastStage
529        % Used f7 for unperturbed solution to interpolate.
530        % Now reset f7 to move along constraint.
531        f(:,nstages) = callODEFunctionNSM( ...
532            ODEFunction,tnew,ynew,options,varargin{:});
533        nfevals = nfevals + 1;
534    end
535    f(:,1) = f(:,nstages); % Already have f(tnew,ynew)
536end
537eml_invariant(nargout <= 3 || (HasEvents && nargout <= 6), ...
538    'MATLAB:maxlhs');
539varargout{1} = tout(1:nout).';
540varargout{2} = yout(:,1:nout).';
541if HasEvents
542    varargout{3} = teout.';
543    varargout{4} = yeout.';
544    varargout{5} = double(ieout.');
545    varargout{6} = double([nsteps;nfailed;nfevals;0;0;0]);
546else
547    varargout{3} = double([nsteps;nfailed;nfevals;0;0;0]);
548    % Spurious definitions required by the compiler.
549    varargout{4} = [];
550    varargout{5} = [];
551    varargout{6} = [];
552end
553if HasOutputFunction
554    OutputFunction(realnull,datanull,'done',varargin{:});
555end
556
557%--------------------------------------------------------------------------
558
559function p = ismonotonic(x)
560coder.internal.prefer_const(x);
561p = true;
562if x(2) > x(1)
563    for k = 2:numel(x)-1
564        if ~(x(k) < x(k+1))
565            p = false;
566            break
567        end
568    end
569else
570    % Start loop off at 1 to catch x(2) == x(1).
571    for k = 1:numel(x)-1
572        if ~(x(k) > x(k+1))
573            p = false;
574            break
575        end
576    end
577end
578
579%--------------------------------------------------------------------------
580
581function p = ispositive(x)
582coder.internal.prefer_const(x);
583p = true;
584for k = 1:numel(x)
585    if ~(x(k) > 0)
586        p = false;
587        break
588    end
589end
590
591%--------------------------------------------------------------------------
592
593function r = normNegative(y,errwt,idxNonNegative,Euclidean)
594% Returns the norm of the negative elements of y(idxNonNegative)./W, where
595% W = errwt, if errwt is a scalar, otherwise, W = errwt(idxNonNegative).
596% Set Euclidean to true to use the 2-norm, otherwise the inf-norm is used.
597coder.inline('always');
598coder.internal.prefer_const(idxNonNegative,Euclidean);
599n = coder.internal.indexInt(numel(idxNonNegative));
600ZERO = zeros('like',real(y));
601if n < 1
602    r = ZERO;
603elseif Euclidean % 2-norm
604    absynn = -real(y(idxNonNegative));
605    for k = 1:n
606        if absynn(k) < 0
607            absynn(k) = ZERO;
608        else
609            absynn(k) = absynn(k)/eml_scalexp_subsref(errwt,idxNonNegative(k));
610        end
611    end
612    r = norm(absynn);
613else % inf-norm
614    r = ZERO;
615    for k = 1:n
616        yk = -real(y(idxNonNegative(k)));
617        if yk > 0
618            yk = yk/eml_scalexp_subsref(errwt,idxNonNegative(k));
619            if yk > r
620                r = yk;
621            end
622        end
623    end
624end
625
626%--------------------------------------------------------------------------
627
628function y = max3(a,b,c)
629% y = max(max(a,b),c);
630% This function requires that c is not NaN.
631y = coder.nullcopy(zeros('like',real(a)));
632if a > b || isnan(b)
633    if a > c
634        y(1) = a;
635    else
636        y(1) = c;
637    end
638elseif b > c
639    y(1) = b;
640else
641    y(1) = c;
642end
643
644%--------------------------------------------------------------------------
645
646function y = maxAbsThresh(a,b,thresh)
647% y = max(abs(a),abs(b),thresh)
648% Inputs a and b must be the same size. Thresh can be scalar or the same
649% size as a and b, but thresh must not be NaN.
650n = coder.internal.indexInt(numel(a));
651y = zeros(size(a),'like',a);
652for k = 1:n
653    y(k) = max3(abs(a(k)),abs(b(k)), ...
654        eml_scalexp_subsref(thresh,k));
655end
656
657%--------------------------------------------------------------------------
658
659function s = make_error_id(solver_name,id)
660coder.inline('always');
661coder.internal.prefer_const(id,solver_name);
662s = coder.const(['MATLAB:',solver_name,':',id]);
663
664%--------------------------------------------------------------------------
665
666function y = REQUESTEDPOINTS
667coder.inline('always');
668y = uint8(0);
669
670function y = SOLVERSTEPS
671coder.inline('always');
672y = uint8(1);
673
674function y = REFINEDSTEPS
675coder.inline('always');
676y = uint8(2);
677
678%--------------------------------------------------------------------------
679
680function next = NextNext(next,ntspan,tspan,tdir,tnew)
681coder.inline('always');
682while next <= ntspan && tdir*(tnew - tspan(next)) >= 0
683    next = next + 1;
684end
685
686%--------------------------------------------------------------------------
687
688function idxNonNegative = getIdxNonNegative(options,y0)
689idx = odeget(options,'NonNegative',[],'fast');
690if isempty(idx)
691    idxNonNegative = zeros(1,0,coder.internal.indexIntClass);
692else
693    neq = coder.internal.indexInt(numel(y0));
694    eml_invariant(isValidIdxVector(idx,neq), ...
695        'MATLAB:odenonnegative:NonNegativeIndicesInvalid');
696    idxNonNegative = coder.internal.indexInt(idx);
697    eml_invariant(isNonNegativeIC(y0,idxNonNegative), ...
698        'MATLAB:odenonnegative:NonNegativeViolatedAtT0');
699end
700
701%--------------------------------------------------------------------------
702
703function p = isValidIdxVector(idx,neq)
704% p = all(idx >= 1 & idx <= neq & idx == floor(idx))
705p = true;
706for k = 1:numel(idx)
707    if ~(idx(k) >= 1 && idx(k) <= neq && idx(k) == floor(idx(k)))
708        p = false;
709        break
710    end
711end
712
713%--------------------------------------------------------------------------
714
715function p = isNonNegativeIC(y0,idxNonNegative)
716% p = all(y0(idxNonNegative) >= 0)
717p = true;
718for k = 1:numel(idxNonNegative)
719    if ~(y0(idxNonNegative(k)) >= 0)
720        p = false;
721        break
722    end
723end
724
725%--------------------------------------------------------------------------
726