MATLAB: Ode23t does not return values for all t in tspan

MATLABode solvertspan

I have a system of ODEs which I am trying to solve numerically as follows:
[t,u]=ode23t(@(ta,AB1234)MyFunc(ta, AB1234,params), [t0 T], [0;0;0;0],options)
where T=0.5. My problem is that the output that I receive does not cover the whole tspan=[t0 T]. For example, I get t values ranging from 0 to 10^-6 rather than [0 0.5]. Does any one know why is that the case? I have already tried other stiff/non-stiff ode solvers. In case of non stiff solver, such as ode45, t contains all values in [0 0.5] but then u is NaN. I have also tried options to reduce 'RelTol' and 'AbsTol', but it didn't help.
clc
clear
warning off
% Initialization
params.v1bar = 0.0084;
params.kappa1 = 9.7196;
params.sigma1 = 0.3924;
params.mu1 = 13.4143;
params.v2bar = 0.0391;
params.kappa2 = 0.1680;
params.sigma2 = 0.1078;
params.mu3 = 0.9238;
params.kappa3 = 0.5967;
params.rho3 = 0.0005;
params.c0neg = 0;
params.c0pos = 1.5713;
params.c1neg = 25.3536;
params.c1pos = 92.4094;
params.c2neg = 0.8802;
params.c2pos = 72.5628;
params.c3neg = 41.4017;
params.c3pos = 0;
params.lambdaneg = 18.7455;
params.lambdapos = 58.2399;
params.eta = 0;
params.delta = 0.01; 
params.r = 0.02;
rho2 = -1;
rho1 =-1;
Y0 =0;
V10 = 0.05;
V20 = 0.03;
V30 = 0.001;
T=0.5;
nu=[0:1:400];
% Calculate the logChar_min_uy
logChar_min_uy_1= logChar_min_uy(1-i*nu,params, T, rho1, rho2);
% Calculate the system of ODEs
function f = logChar_min_uy(u,params, T, rho1, rho2)
t0=0;
options = odeset('AbsTol',1e-9,'Refine',6);
for ii=1:length(u)
[Tnumer, alpha_beta1234Numeric] = ode23t(@(ta,AB1234)alpha_beta1234(ta, AB1234,u(ii),params, rho1, rho2), [t0 T], [0;0;0;0],options);
f(ii,:) = alpha_beta1234Numeric(end,:);
end
end
% System of ODEs
function Dalphabeta_Dtau = alpha_beta1234(tau, AB1234,u,params, rho1, rho2)
alpha = AB1234(1);
beta1 = AB1234(2);
beta2 = AB1234(3);
beta3 = AB1234(4);
% Calculate the integrals
Thetanc_u00 = integral(@(z) params.lambdaneg*exp(z*(u+params.lambdaneg)), -inf,0);
Thetanc_ubeta1beta3 = integral(@(z) params.lambdaneg*exp(z*(u+params.lambdaneg)+beta1*params.mu1*z.^2+beta3*(1-params.rho3)*params.mu3*z.^2), -inf,0);
Thetap_u = integral(@(z) params.lambdapos*exp(z*(u-params.lambdapos)), 0,inf);
Thetani_beta3 =  integral(@(z) params.lambdaneg*exp(z*params.lambdaneg+beta3*params.rho3*params.mu3*z.^2), -inf,0);
% ODEs
dalpha_dtau = u*(params.r- params.delta- params.c0neg*(Thetanc_u00-1) - params.c0pos*(Thetap_u-1))+ beta1*params.kappa1*params.v1bar                                    + beta2*params.kappa2*params.v2bar + params.c0neg*(Thetanc_ubeta1beta3-1) + params.c0neg*(Thetani_beta3-1) + params.c0pos*(Thetap_u-1); 
dbeta1_dtau = u*(-0.5                  - params.c1neg*(Thetanc_u00-1) - params.c1pos*(Thetap_u-1))- beta1*params.kappa1 + 0.5 *u^2            +0.5* params.sigma1^2*beta1^2 + beta1 * u*params.sigma1* rho1+ params.c1neg*(Thetanc_ubeta1beta3-1) + params.c1neg*(Thetani_beta3-1) + params.c1pos*(Thetap_u-1); 
dbeta2_dtau = u*(-0.5                  - params.c2neg*(Thetanc_u00-1) - params.c2pos*(Thetap_u-1))- beta2*params.kappa2 + 0.5 *u^2            +0.5* params.sigma2^2*beta2^2 + beta2 * u*params.sigma2* rho2+ params.c2neg*(Thetanc_ubeta1beta3-1) + params.c2neg*(Thetani_beta3-1) + params.c2pos*(Thetap_u-1);
dbeta3_dtau = u*(-0.5*params.eta^2     - params.c3neg*(Thetanc_u00-1) - params.c3pos*(Thetap_u-1))- beta3*params.kappa3 + 0.5 *u^2*params.eta^2                                                            + params.c3neg*(Thetanc_ubeta1beta3-1) + params.c3neg*(Thetani_beta3-1) + params.c3pos*(Thetap_u-1);
Dalphabeta_Dtau = [dalpha_dtau;dbeta1_dtau;dbeta2_dtau;dbeta3_dtau];
end

Best Answer

You have
Thetanc_ubeta1beta3 = integral(@(z) params.lambdaneg*exp(z*(u+params.lambdaneg)+beta1*params.mu1*z.^2+beta3*(1-params.rho3)*params.mu3*z.^2), -inf,0);
Notice there is exp() of an expression that includes beta1*z^2 .
At some point, your boundary conditions beta1 = AB1234(2) becomes complex valued, so you have the integral of a complex value times z^2 from -inf to 0. The theoretical value for that integral is not defined. If you ask MATLAB's Symbolic Toolbox for the integral, the result will be in terms of a limit at infinity that is not defined. Maple (a different software package) just plain says that the integral is not defined.
When you ask for numeric integration of this undefined integral, rather than returning empty or generating an error, MATLAB returns complex nan, nan+1i*nan . That then polutes the rest of your calculations.
Related Question