function [tout,yout] = ode45(ode,tspan,y0,options,varargin)
%ODE45  Solve non-stiff differential equations, medium order method.
%   [TOUT,YOUT] = ODE45(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL]
%   integrates the system of differential equations y' = f(t,y) from time
%   T0 to TFINAL with initial conditions Y0. ODEFUN is a function handle.
%   For a scalar T and a vector Y, ODEFUN(T,Y) must return a column vector
%   corresponding to f(t,y). Each row in the solution array YOUT
%   corresponds to a time returned in the column vector TOUT. To obtain
%   solutions at specific times T0,T1,...,TFINAL (all increasing or all
%   decreasing), use TSPAN = [T0 T1 ... TFINAL]. At least one between Y0
%   and TSPAN must be an unformatted dlarray object. If Y0 or TSPAN is a
%   dlarray object, then also YOUT and TOUT is a dlarray, respectively.
%   
%   [TOUT,YOUT] = ODE45(ODEFUN,TSPAN,Y0,OPTIONS) solves as above with default
%   integration properties replaced by values in OPTIONS, an argument created
%   with the ODESET function. See ODESET for details. Commonly used options 
%   are scalar relative error tolerance 'RelTol' (1e-3 by default) and vector
%   of absolute error tolerances 'AbsTol' (all components 1e-6 by default).
%
%   For dlarray input, the ode45 function does not support the OutputFcn,
%   Mass, NonNegative, and Events options.
% 
%   Example    
%         [t,y]=ode45(@vdp1,[0 20],dlarray([2 0]));   
%     solves the system y' = vdp1(t,y), using the default relative error
%     tolerance 1e-3 and the default absolute tolerance of 1e-6 for each
%     component, and plots the first component of the solution. 
%   
%   Class support for inputs TSPAN, Y0, and the result of ODEFUN(T,Y):
%     float: double, single, dlarray
%
%   See also DLARRAY, DLODE45.

%   ODE45 is an implementation of the explicit Runge-Kutta (4,5) pair of
%   Dormand and Prince called variously RK5(4)7FM, DOPRI5, DP(4,5) and DP54.
%   It uses a "free" interpolant of order 4 communicated privately by
%   Dormand and Prince.  Local extrapolation is done.

%   Details are to be found in The MATLAB ODE Suite, L. F. Shampine and
%   M. W. Reichelt, SIAM Journal on Scientific Computing, 18-1, 1997.

%   Mark W. Reichelt and Lawrence F. Shampine, 6-14-94
%   Copyright 1984-2021 The MathWorks, Inc.

arguments
   ode
   tspan {mustBeVector}
   y0 {mustBeA(y0,'dlarray'),mustBeReal}
   options = odeset()
end
arguments(Repeating)
   varargin
end

isTspanDlarray = isdlarray(tspan);

if nargout ~= 2
    error(message("deep:dlarray:Ode45TwoOutputSyntax"));
end

if ~isempty(varargin)
    error(message("deep:dlarray:Ode45AdditionalInput"));
end

if (isdlarray(y0) && ~isempty(dims(y0))) || (isTspanDlarray && ~isempty(dims(tspan)) )
    error(message("deep:dlarray:Ode45format"));
end

if isgpuarray(y0) || isgpuarray(tspan)
    error(message("deep:dlarray:Ode45gpuArray"));
end
    
iValidateIsDefaultOption(options,'Mass');
iValidateIsDefaultOption(options,'Events');
iValidateIsDefaultOption(options,'NonNegative');
iValidateIsDefaultOption(options,'OutputFcn');

solver = 'ode45';

if isempty(tspan) || isempty(y0) 
    error(message('MATLAB:odearguments:TspanOrY0NotSupplied', solver));
end

if length(tspan) < 2
    error(message('MATLAB:odearguments:SizeTspan', solver));
end

isY0Recorded    = ~isempty(getRecorded(y0));
isTspanRecorded = isTspanDlarray && ~isempty(getRecorded(tspan));

if isY0Recorded || isTspanRecorded
    [y0,fd] = deep.internal.dlarray.toInternal(y0);
    % Extract for initial condition data.
    [fd, y0data] = extractData(fd);
    
    if isTspanDlarray
        [tspan,fd_tspan] = deep.internal.dlarray.toInternal(tspan);
        [~,tspan_data] = extractData(fd_tspan);
    else
        tspan_data = tspan;
    end

    % Record direct identity operation.
    [tspan_recording,y0_recording] = internal_dlode_direct_identity(tspan_data,y0data,"ode45");

    if isY0Recorded
        fd.Layout = deep.internal.format.DataLayout('',[]);
        fd = setDataAndLayout(fd, matlab.lang.internal.move(y0_recording), fd.Layout);
        y0 = deep.internal.dlarray.fromInternal(fd,y0);
    end

    if isTspanDlarray
        fd_tspan.Layout = deep.internal.format.DataLayout('',[]);
        fd_tspan = setDataAndLayout(fd_tspan, matlab.lang.internal.move(tspan_recording), fd_tspan.Layout);
        tspan = deep.internal.dlarray.fromInternal(fd_tspan,tspan);
    end
end

if isdlarray(tspan)
  t0 = parenReference(tspan,1);
  t1 = parenReference(tspan,2);
  tfinal = parenReference(tspan,numel(tspan));
else
  t0 = tspan(1);
  t1 = tspan(2);
  tfinal = tspan(end);
end
if t0 == tfinal
    error(message('MATLAB:odearguments:TspanEndpointsNotDistinct'));
end

y0 = reshape(y0,[],1);
neq = length(y0);

FcnHandlesUsed = isa(ode,'function_handle');
if ~FcnHandlesUsed
    % ode-file used. This is not an accepted syntax for dlarray overload of
    % ode45.
    error(message("deep:dlarray:Ode45odeFile"));
end

if nargin(ode) ~= 2
    error(message("deep:dlarray:OdeOdefunTwoInputs"));
end

f0 = feval(ode,t0,y0);
[m,n] = size(f0);
if n > 1
    error(message('MATLAB:odearguments:FoMustReturnCol', func2str( ode )));
elseif m ~= neq
    error(message('MATLAB:odearguments:SizeIC', func2str( ode ), m, neq, func2str( ode )));
end

% Determine the dominant data type
classT0 = underlyingType(t0);
classY0 = underlyingType(y0);
classF0 = underlyingType(f0);

dataType = superiorfloat(iExtractIfDlarray(tspan),iExtractIfDlarray(y0),iExtractIfDlarray(f0));

if ~(strcmp(classT0,dataType) && strcmp(classY0,dataType) && strcmp(classF0,dataType))
    input1 = '''t0'', ''y0'''; 
    input2 = '''f(t0,y0)''';
    warning(message('MATLAB:odearguments:InconsistentDataType',input1,input2,solver));
end        

t0 = cast(t0, dataType);
y0 = cast(y0, dataType);
f0 = cast(f0, dataType);

% Call function implementing the numerical integration loop.
[tout,yout] = nnet.internal.cnnhost.ode45(ode,tspan,y0,options,t0,t1,tfinal,f0);
end

function iValidateIsDefaultOption(options, optionName)
if ~isequal( options.(optionName), [] )
    error(message("deep:dlarray:Ode45NotDefault",optionName));
end
end

function x = iExtractIfDlarray(x)
if isdlarray(x)
    x = extractdata(x);
end
end
